Metagenomic analysis of Limnohabitans population in Lake Michigan
Ruben Props
Today
# Read data
mean_coverage <- read.table("./anvio_output/rebin/mean_coverage_selected_final.tsv", header = TRUE)
std_coverage <- read.table("./anvio_output/rebin/std_coverage_selected_final.tsv", header = TRUE)
bin_size <- read.table("./anvio_output/rebin/general_bins_summary_selected_final.tsv", header = TRUE)[, c(2,5)]
total_reads <- read.table("./anvio_output/sample_reads.tsv", header = TRUE)
read_length <- 150
# From wide to long format
mean_coverage_long <- tidyr::gather(mean_coverage, Sample_ID, coverage,
Fa13_BD_MLB_DN:Su13_BD_MM15_SN_C, factor_key=TRUE)
mean_coverage_long[,2] <- gsub(mean_coverage_long[,2], pattern = "_C",
replacement = "")
std_coverage_long <- tidyr::gather(std_coverage, Sample_ID, std_coverage,
Fa13_BD_MLB_DN:Su13_BD_MM15_SN_C,
factor_key=TRUE)
std_coverage_long[,2] <- gsub(std_coverage_long[,2], pattern = "_C",
replacement = "")
coverage_data <- data.frame(mean_coverage_long,
std_coverage = std_coverage_long[,3])
# Read and add metadata
meta <- read.csv2("metadata.csv")
meta$Sample_ID <- gsub(meta$Sample_ID, pattern = ".", replacement = "_", fixed = TRUE)
data_total <- left_join(coverage_data, total_reads, by = c("Sample_ID" = "sample"))
data_total <- left_join(data_total, bin_size, by = "bins")
data_total <- left_join(data_total, meta, by = "Sample_ID")
# Calculate relative abundance of the bins
data_total$mean_rel_abundance <- 100*(data_total$coverage*data_total$bin_size)/(read_length*data_total$Total_reads)
data_total$upper_rel_abundance <- 100*((data_total$coverage+data_total$std_coverage)*data_total$bin_size)/(read_length*data_total$Total_reads)
data_total$lower_rel_abundance <- 100*((data_total$coverage-data_total$std_coverage)*data_total$bin_size)/(read_length*data_total$Total_reads)1a. Phylogeny
Phylogenomic tree
16S rRNA gene phylogenetic tree
Annotated 16S rRNA gene phylogenetic tree
Microdiversity in MAG8
At each iteration, SSU sequences were merged into one sequence if the identity of non-gapped positions in a global alignment was greater than 97%. A single SSU sequence (and its prior probability) was divided into two sequences if the second most probable base in more than 4% of all positions had a probability greater than 10%. In this way, sequences that evolved over iterations to be the same were merged, and sequences with evidence from the reads for multiple OTUs were duplicated and allowed to evolve as separate OTUs in future iterations. Reference
The workflow for postprocessing the EMIRGE reconstructed 16S rRNA gene sequences:
# Import seq-table that was generated after EMIRGE post-processing
EMIRGE_OTU <- read.table("./Tree/EMIRGE/read.info.txt",
stringsAsFactors = FALSE)
colnames(EMIRGE_OTU) <- c("Sample", "Sequence", "Prior", "Length", "Prior_norm")
EMIRGE_names <- read.table("./Tree/EMIRGE/EMIRGE_Limno.names")
EMIRGE_OTU_97 <- read.table("./Tree/EMIRGE/final.otu.emirge.pick.an.0.03.rep.names",
stringsAsFactors = FALSE)
EMIRGE_OTU_99 <- read.table("./Tree/EMIRGE/final.otu.emirge.pick.an.0.01.rep.names",
stringsAsFactors = FALSE)
# Filter out putative Limnohabitans OTUs
EMIRGE_OTU <- EMIRGE_OTU[EMIRGE_OTU$Sequence %in% EMIRGE_names$V1, ]
# Formate sample names
EMIRGE_OTU$Sample <- gsub(".A","", EMIRGE_OTU$Sample, fixed = TRUE)
EMIRGE_OTU$Sample <- gsub(".C","", EMIRGE_OTU$Sample, fixed = TRUE)
# Add OTU cluster results as new columns (0.97 level and 0.99 level)
EMIRGE_OTU$OTU_97 <- EMIRGE_OTU$Sequence
for(i in 1:nrow(EMIRGE_OTU_97)){
EMIRGE_OTU$OTU_97[EMIRGE_OTU$Sequence %in%
do.call(rbind, strsplit(EMIRGE_OTU_97$V2[i], ","))] <-
EMIRGE_OTU_97$V3[i]
}
EMIRGE_OTU$OTU_99 <- EMIRGE_OTU$Sequence
for(i in 1:nrow(EMIRGE_OTU_99)){
EMIRGE_OTU$OTU_99[EMIRGE_OTU$Sequence %in%
do.call(rbind, strsplit(EMIRGE_OTU_99$V2[i], ","))] <-
EMIRGE_OTU_99$V3[i]
}
# Calculate new abundance for the relative abundance of each OTU (0.97 and 0.99 level)
EMIRGE_OTU <- EMIRGE_OTU %>% dplyr::group_by(OTU_97, Sample) %>%
dplyr::mutate(Prior_norm_OTU97 = sum(Prior_norm))
EMIRGE_OTU <- EMIRGE_OTU %>% dplyr::group_by(OTU_99, Sample) %>%
dplyr::mutate(Prior_norm_OTU99 = sum(Prior_norm))
# Add metadata variables as new column
meta_em <- meta[-1]; meta_em$Sample_ID <- gsub("_",".", meta_em$Sample_ID)
meta_em <- meta_em %>% dplyr::distinct()
EMIRGE_OTU <- dplyr::left_join(EMIRGE_OTU, meta_em, by = c("Sample" = "Sample_ID"))
# to change# Plot
p_emirge_abund_otu97 <- ggplot(EMIRGE_OTU, aes(x = Sample, y = 100*Prior_norm, fill = OTU_97))+
theme_bw()+
geom_bar(stat = "identity", size = 1, color = "black", alpha = 0.7)+
scale_fill_brewer(palette = "Accent")+
theme(axis.text=element_text(size=14), axis.title=element_text(size=20),
axis.text.x = element_text(size = 14, angle = 45, hjust = 1),
title=element_text(size=20), legend.text=element_text(size=12),
legend.background = element_rect(fill="transparent"),
strip.text=element_text(size=14), legend.position = "right",
strip.background = element_rect(fill = adjustcolor("gray", 0.15)))+
guides(fill=guide_legend(nrow = 3))+
facet_grid(.~Site, scales = "free_x")+
ylab("Length-normalized prior (%)")+
xlab("")+
ggtitle("97% clustered OTUs")
# coord_trans(x = "atanh", y = "atanh")
# Plot
p_emirge_abund_otu99 <- ggplot(EMIRGE_OTU, aes(x = Sample, y = 100*Prior_norm, fill = OTU_99))+
theme_bw()+
geom_bar(stat = "identity", size = 1, color = "black", alpha = 0.7)+
# geom_point(size = 4, shape = 21, color = "black", alpha = 0.7)+
scale_fill_brewer(palette = "Accent")+
theme(axis.text=element_text(size=14), axis.title=element_text(size=20),
axis.text.x = element_text(size = 14, angle = 45, hjust = 1),
title=element_text(size=20), legend.text=element_text(size=12),
legend.background = element_rect(fill="transparent"),
strip.text=element_text(size=14), legend.position = "right",
strip.background = element_rect(fill = adjustcolor("gray", 0.15)))+
guides(fill=guide_legend(nrow = 3))+
facet_grid(.~Site, scales = "free_x")+
ylab("Length-normalized prior (%)")+
xlab("")+
ggtitle("99% clustered OTUs")
# coord_trans(x = "atanh", y = "atanh")
# Plot
p_emirge_abund_seq <- ggplot(EMIRGE_OTU, aes(x = Sample, y = 100*Prior_norm, fill = OTU_97))+
theme_bw()+
# geom_point(size = 4, shape = 21, color = "black", alpha = 0.7)+
geom_bar(stat = "identity", size = 1, color = "black", alpha = 0.7)+
scale_fill_brewer(palette = "Accent")+
theme(axis.text=element_text(size=14), axis.title=element_text(size=20),
axis.text.x = element_text(size = 14, angle = 45, hjust = 1),
title=element_text(size=20), legend.text=element_text(size=12),
legend.background = element_rect(fill="transparent"),
strip.text=element_text(size=14), legend.position = "right",
strip.background = element_rect(fill = adjustcolor("gray", 0.15)))+
guides(fill=guide_legend(nrow = 3))+
facet_grid(.~Site, scales = "free_x")+
ylab("Length-normalized prior (%)")+
xlab("")+
ggtitle("Unique sequences")
# coord_trans(x = "atanh", y = "atanh")
grid.arrange(p_emirge_abund_otu97, p_emirge_abund_otu99, p_emirge_abund_seq ,nrow = 3)
# Import blast results
blast_emirge <- read.csv("./Tree/EMIRGE/blast_table.csv", header = FALSE)[, 1:3]
colnames(blast_emirge) <- c("seq1", "seq2", "pid")
# Annotate based on 97% clusters
EMIRGE_OTU_annot <- EMIRGE_OTU %>% dplyr::select(Sequence, OTU_97) %>%
distinct() %>%
left_join(., EMIRGE_OTU_97[, c("V3", "V4")], by = c("OTU_97" = "V3"))## Adding missing grouping variables: `OTU_99`, `Sample`# label sequence 1
blast_emirge_annot <- left_join(blast_emirge, EMIRGE_OTU_annot,
by = c("seq1" = "Sequence")) %>%
dplyr::filter(V4 != "not_limno")
colnames(blast_emirge_annot)[7] <- "seq1_label"
# label sequence 2
blast_emirge_annot <- left_join(blast_emirge_annot, EMIRGE_OTU_annot[, c(3,5)],
by = c("seq2" = "Sequence")) %>%
dplyr::filter(V4 != "not_limno")
colnames(blast_emirge_annot)[8] <- "seq2_label"
# Create interaction effect
blast_emirge_annot$label_interaction <- interaction(blast_emirge_annot$seq1_label,
blast_emirge_annot$seq2_label,
sep= "-")
# Only visualize within-group blast identities
subset_seq_labels <- c("LimB-LimB", "Lhab-A4-Lhab-A4")
blast_emirge_annot %>%
dplyr::filter(label_interaction %in% subset_seq_labels &
seq1 != seq2) %>%
ggplot(aes(x = label_interaction, y = pid, fill = label_interaction))+
geom_jitter(size = 2.5,
color = "black", shape = 21, width = 0.1, alpha = 0.5)+
geom_violin(alpha = 0.4, adjust = 1, draw_quantiles = TRUE, width = 0.3)+
stat_summary(fun.data=mean_sdl, fun.args = list(mult = 1),
geom="pointrange", color="black")+
ylim(95,100)+
scale_fill_manual("", values = rev(c("#e31a1c","#fdbf6f")))+
theme_bw()+
theme(axis.text=element_text(size=14), axis.title=element_text(size=20),
axis.text.x=element_text(size=14),
title=element_text(size=16), legend.text=element_text(size=14),
legend.background = element_rect(fill="transparent"),
strip.text.x=element_text(size=18))+
guides(fill = FALSE)+
labs(x="")
blast_emirge_annot %>% dplyr::filter(label_interaction %in% subset_seq_labels &
seq1 != seq2) %>%
group_by(label_interaction) %>%
summarize(mean(pid), sd(pid))## # A tibble: 2 x 3
## label_interaction `mean(pid)` `sd(pid)`
## <fct> <dbl> <dbl>
## 1 Lhab-A4-Lhab-A4 98.8 0.826
## 2 LimB-LimB 98.8 0.850# Only look at the outlying sequence of LimB
# 10|JN626652.1.1352
blast_emirge_annot$label_interaction <- interaction(blast_emirge_annot$seq1_label,
blast_emirge_annot$seq2_label,
sep= "-")
# Only visualize within-group blast identities
subset_seq_labels <- c("LimB-LimB", "Lhab-A4-Lhab-A4")
blast_emirge_annot %>% dplyr::filter(label_interaction %in% subset_seq_labels &
seq1 == "10|JN626652.1.1352") %>%
summarize(mean(pid), sd(pid))## mean(pid) sd(pid)
## 1 97.96535 0.6247723ANI analysis
# read file
data.ANI <- read.table("./ANI/ANIb_percentage_identity.tab")
# Replace genome names by better annotated names
map_ani <- read.delim("./ANI/pyani_ref_names.tsv", stringsAsFactors = FALSE)
for(i in 1:nrow(map_ani)){
colnames(data.ANI)[colnames(data.ANI) %in% map_ani$ani_file[i]] <- map_ani$ref_name[i]
rownames(data.ANI)[rownames(data.ANI) %in% map_ani$ani_file[i]] <- map_ani$ref_name[i]
}
# Remove pnec and root genomes from dataframe
data.ANI <- data.ANI[grep("MAG|Lim", rownames(data.ANI)), grep("MAG|Lim", rownames(data.ANI))]
# Order y-axis according to phylogenetic tree order
ord_list_bin <- c(
"MAG5.SP-M110-DD", "MAG2.FA-MLB-SN",
"MAG3.FA-MLB-SN", "MAG4.FA-M110-DN",
"MAG1.FA-MLB-DN", "MAG10.SU-M15-SN",
"MAG6.SP-M15-SD","MAG8.SU-M110-DCMD",
"Lim. sp. Rim11", "Lim. sp. 103DPR2",
"MAG7.SU-MLB-SD", "Lim. sp. Rim47",
"MAG9.SU-M15-SN", "Lim. sp. Rim28",
"Lim. sp. 63ED37-2","Lim. sp. 2KL-27",
"Lim. sp. 2KL-3", "Lim. sp. DM1",
"Lim. sp. II-D5"
)
# order rows and columns
data.ANI <- data.ANI[, order(match(colnames(data.ANI), ord_list_bin))]
data.ANI <- data.ANI[order(match(rownames(data.ANI), ord_list_bin)), ]
data.ANI <- get_upper_tri(data.ANI)
# melt to dataframe
df_pyani <- melt(as.matrix(data.ANI), na.rm = TRUE) # reshape into dataframe
names(df_pyani)[c(1:3)] <- c("bin1", "bin2", "ANI")
df_pyani[, 1] <- as.character(df_pyani[, 1]); df_pyani[, 2] <- as.character(df_pyani[, 2])
df_pyani$bin2 <- factor(df_pyani$bin2, levels = ord_list_bin)
df_pyani$bin1 <- factor(df_pyani$bin1, levels = rev(ord_list_bin))
# make plot
p_ani <- ggplot(data = df_pyani,
aes(x = bin2, y = bin1, fill = 100*ANI))+
geom_raster()+
scale_fill_distiller(palette = "RdBu", name = "Average\nNucleotide\nIdentity (%)\n",
limits = 100*c(0.75,1.0), oob=squish) +
geom_text(aes(label = round(100*ANI, 0)), size = 4.5)+
xlab("")+
ylab("")+
scale_x_discrete(position = "top") +
theme(axis.title=element_text(size=14), strip.text.x=element_text(size=14),
legend.title=element_text(size=14), legend.text=element_text(size=14),
axis.text.y = element_blank(),
axis.text.x = element_text(size=14, angle = 55, hjust = 0),
title=element_text(size=20),
plot.margin = unit(c(1.1,1.1,1.1,1.1), "cm"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_rect(fill = "transparent",colour = NA),
plot.background = element_rect(fill = "transparent",colour = NA)
)
# png("ANI_heatmap.png", height= 9, width = 9, res = 500, units = "in", bg = "transparent")
print(p_ani)
# dev.off()1b. Investigate predicted generation time (MGT)
MGT_df <- read.table("./Growthpred/GP_results.tsv", header = TRUE,
stringsAsFactors = FALSE, sep = "\t")
# Order genome_ids according to the phylogenetic clustering
MGT_df$Genome_ID <- factor(MGT_df$Genome_ID, levels = ord_list_bin)
MGT_df$Lineage <- factor(MGT_df$Lineage, levels = c("LimDEA_1", "LimDEA_2",
"LimB", "LimC"))
# Make barplot with st.dev to visualize MGT and optimal temperature
selected_points <- data.frame(Genome_ID = MGT_df$Genome_ID,
ypos = c(rep(0.25, 10), rep(NA,9)))
p_MGT <- ggplot(MGT_df, aes(x = Genome_ID, y = log(2)/MGT, fill = Lineage, group = Genome_ID))+
theme_bw()+
geom_bar(alpha = 1, stat = "identity", color = "black",
position = position_dodge(width = 1), width = 0.7)+
scale_fill_manual(expression("Growth rate - "~"h"^-1),
values = c(rgb(red=t(col2rgb("#deebf7ff")), maxColorValue = 255),
rgb(red=t(col2rgb("#c6dbefff")), maxColorValue = 255),
rgb(red=t(col2rgb("#9ecae1ff")), maxColorValue = 255),
rgb(red=t(col2rgb("#6baed6ff")), maxColorValue = 255)
))+
theme(axis.text=element_text(size=15, face = "bold"), axis.title=element_text(size=20),
title=element_text(size=20), legend.text=element_text(size=14),
legend.background = element_rect(fill="transparent"),
# axis.text.x = element_text(angle = 65, hjust = 1),
strip.text.x=element_text(size = 18),
legend.position="bottom",
axis.text.x=element_text(size = 13, angle =45, hjust= 1),
axis.title.x=element_blank(),
plot.title = element_text(hjust = 0, size=18))+
guides(fill=FALSE)+
geom_errorbar(aes(ymin = log(2)/(MGT - sd.MGT), ymax = log(2)/(MGT + sd.MGT)), width = 0.15,
position = position_dodge(width = 1))+
ylab("")+
xlab("")+
ggtitle(expression("Growth rate - h"^-1))+
geom_point(data = selected_points, aes(x = Genome_ID, y = ypos),
shape = 25, fill = "black", col = "black", size = 3)
# svg(filename = "MGT_figure.svg", width = 9.5, height = 4)
print(p_MGT)
# dev.off()
# Make the same plot for optimal growth temperature
selected_points <- data.frame(Genome_ID = MGT_df$Genome_ID,
ypos = c(rep(30, 10), rep(NA,9)))
p_Topt <- ggplot(MGT_df, aes(x = Genome_ID, y = Topt, fill = Lineage, group = Genome_ID))+
theme_bw()+
geom_bar(alpha = 0.4, stat = "identity", color = "black",
position = position_dodge(width = 1), width = 0.7)+
scale_fill_manual("Optimal growth temperature (°C)",
values = c("#deebf7ff", "#c6dbefff","#9ecae1ff",
"#6baed6ff"))+
theme(axis.text=element_text(size=15, face = "bold"), axis.title=element_text(size=20),
title=element_text(size=20), legend.text=element_text(size=14),
legend.background = element_rect(fill="transparent"),
# axis.text.x = element_text(angle = 65, hjust = 1),
strip.text.x=element_text(size = 18),
legend.position="bottom",
axis.text.x=element_text(size = 13, angle =45, hjust= 1),
axis.title.x=element_blank(),
plot.title = element_text(hjust = 0, size=18))+
guides(fill=FALSE)+
ylab("")+
xlab("")+
ggtitle("Optimal growth temperature (°C)")+
geom_point(data = selected_points, aes(x = Genome_ID, y = ypos),
shape = 25, fill = "black", col = "black", size = 3)
print(p_Topt)
1c. Network analysis based on 16S data
# import data
df_phy <- import_mothur(mothur_shared_file = "./16S/stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.pick.an.unique_list.shared",
mothur_constaxonomy_file = "./16S/stability.trim.contigs.good.unique.good.filter.unique.precluster.pick.pick.an.unique_list.0.03.cons.taxonomy")
# Filter out 2013 samples
df_phy <- prune_samples(grep(pattern = ".", sample_names(df_phy), fixed = TRUE,
value = TRUE), df_phy)
df_phy <- prune_samples(grep(pattern = "cD", sample_names(df_phy), fixed = TRUE,
value = TRUE, invert = TRUE), df_phy)
df_phy <- prune_samples(grep(pattern = "ED", sample_names(df_phy), fixed = TRUE,
value = TRUE, invert = TRUE), df_phy)
# Calculate diversity
diversity_df_phy <- Diversity_16S(df_phy, brea = FALSE, R = 100)## **WARNING** this functions assumes that rows are samples and columns
## are taxa in your phyloseq object, please verify.
## Fri May 25 17:37:00 2018 Calculating diversity for sample 1/41 --- Fa13.BD.MLB.DN.1.renamed
## Fri May 25 17:37:35 2018 Calculating diversity for sample 2/41 --- Fa13.BD.MLB.SN.1.renamed
## Fri May 25 17:38:11 2018 Calculating diversity for sample 3/41 --- Fa13.BD.MM110.DN.1.renamed
## Fri May 25 17:38:44 2018 Calculating diversity for sample 4/41 --- Fa13.BD.MM110.DN.2.renamed
## Fri May 25 17:39:18 2018 Calculating diversity for sample 5/41 --- Fa13.BD.MM110.SD.1.renamed
## Fri May 25 17:39:51 2018 Calculating diversity for sample 6/41 --- Fa13.BD.MM110.SD.2.renamed
## Fri May 25 17:40:27 2018 Calculating diversity for sample 7/41 --- Fa13.BD.MM110.SN.1.renamed
## Fri May 25 17:41:00 2018 Calculating diversity for sample 8/41 --- Fa13.BD.MM110.SN.2.renamed
## Fri May 25 17:41:34 2018 Calculating diversity for sample 9/41 --- Fa13.BD.MM15.DN.1.renamed
## Fri May 25 17:42:06 2018 Calculating diversity for sample 10/41 --- Fa13.BD.MM15.DN.2.renamed
## Fri May 25 17:42:40 2018 Calculating diversity for sample 11/41 --- Fa13.BD.MM15.SD.1.renamed
## Fri May 25 17:43:18 2018 Calculating diversity for sample 12/41 --- Fa13.BD.MM15.SD.2.renamed
## Fri May 25 17:43:51 2018 Calculating diversity for sample 13/41 --- Fa13.BD.MM15.SN.1.renamed
## Fri May 25 17:44:21 2018 Calculating diversity for sample 14/41 --- Fa13.BD.MM15.SN.2.renamed
## Fri May 25 17:44:53 2018 Calculating diversity for sample 15/41 --- Sp13.BD.MLB.SN.1.renamed
## Fri May 25 17:45:29 2018 Calculating diversity for sample 16/41 --- Sp13.BD.MLB.SN.2.renamed
## Fri May 25 17:46:06 2018 Calculating diversity for sample 17/41 --- Sp13.BD.MM110.DD.1.renamed
## Fri May 25 17:46:46 2018 Calculating diversity for sample 18/41 --- Sp13.BD.MM110.SD.1.renamed
## Fri May 25 17:47:25 2018 Calculating diversity for sample 19/41 --- Sp13.BD.MM110.SD.2.renamed
## Fri May 25 17:48:04 2018 Calculating diversity for sample 20/41 --- Sp13.BD.MM110.SN.1.renamed
## Fri May 25 17:48:43 2018 Calculating diversity for sample 21/41 --- Sp13.BD.MM110.SN.2.renamed
## Fri May 25 17:49:27 2018 Calculating diversity for sample 22/41 --- Sp13.BD.MM15.DD.1.renamed
## Fri May 25 17:50:09 2018 Calculating diversity for sample 23/41 --- Sp13.BD.MM15.SD.1.renamed
## Fri May 25 17:50:46 2018 Calculating diversity for sample 24/41 --- Sp13.BD.MM15.SN.1.renamed
## Fri May 25 17:51:24 2018 Calculating diversity for sample 25/41 --- Sp13.BD.MM15.SN.2.renamed
## Fri May 25 17:52:03 2018 Calculating diversity for sample 26/41 --- Su13.BD.MLB.DD.1.renamed
## Fri May 25 17:52:37 2018 Calculating diversity for sample 27/41 --- Su13.BD.MLB.SD.1.renamed
## Fri May 25 17:53:13 2018 Calculating diversity for sample 28/41 --- Su13.BD.MM110.DCMD.1.renamed
## Fri May 25 17:53:53 2018 Calculating diversity for sample 29/41 --- Su13.BD.MM110.DCMD.2.renamed
## Fri May 25 17:54:37 2018 Calculating diversity for sample 30/41 --- Su13.BD.MM110.DN.1.renamed
## Fri May 25 17:55:13 2018 Calculating diversity for sample 31/41 --- Su13.BD.MM110.DN.2.renamed
## Fri May 25 17:55:46 2018 Calculating diversity for sample 32/41 --- Su13.BD.MM110.SD.1.renamed
## Fri May 25 17:56:20 2018 Calculating diversity for sample 33/41 --- Su13.BD.MM110.SD.2.renamed
## Fri May 25 17:56:54 2018 Calculating diversity for sample 34/41 --- Su13.BD.MM110.SN.1.renamed
## Fri May 25 17:57:28 2018 Calculating diversity for sample 35/41 --- Su13.BD.MM110.SN.2.renamed
## Fri May 25 17:58:01 2018 Calculating diversity for sample 36/41 --- Su13.BD.MM15.DN.1.renamed
## Fri May 25 17:58:34 2018 Calculating diversity for sample 37/41 --- Su13.BD.MM15.DN.2.renamed
## Fri May 25 17:59:08 2018 Calculating diversity for sample 38/41 --- Su13.BD.MM15.SD.1.renamed
## Fri May 25 17:59:41 2018 Calculating diversity for sample 39/41 --- Su13.BD.MM15.SD.2.renamed
## Fri May 25 18:00:15 2018 Calculating diversity for sample 40/41 --- Su13.BD.MM15.SN.1.renamed
## Fri May 25 18:00:50 2018 Calculating diversity for sample 41/41 --- Su13.BD.MM15.SN.2.renamed
## Fri May 25 18:01:24 2018 Done with all 41 samples# Perform prevalence filtering
df_phy <- filter_taxa(df_phy, function(x) sum(x > 30) > (0.25*length(x)), TRUE)
# Run spiec-easi
sp_easi <- spiec.easi(df_phy, method='mb', lambda.min.ratio=1e-2,
nlambda=20, icov.select.params=list(rep.num=50))## Normalizing/clr transformation of data with pseudocount ...## Inverse Covariance Estimation with mb ...## Model selection with stars ...## Done!ig.mb <- adj2igraph(sp_easi$refit, vertex.attr = list(name=taxa_names(df_phy)))
vsize <- Biobase::rowMax(clr(t(otu_table(df_phy)), 1))+10
Lineage_rel <- tax_table(df_phy)[,"Rank2"]
Lineage_rel <- factor(Lineage_rel, levels = unique(Lineage_rel))
# OTUs that are Limnohabitans
limno_otus <- taxa_names(df_phy)[tax_table(df_phy)[,"Rank6"] == "betI_A"]
limno_otus <- limno_otus[!is.na(limno_otus)]
# Make Limno label
limno_labs <- c()
limno_labs[vertex.attributes(ig.mb)$name %in% limno_otus] <- "Limnohabitans sp."
limno_labs[is.na(limno_labs)] <- ""
# Plot network
p_16S_network <- plot_network_custom(ig.mb, df_phy, type='taxa',
line_weight = 2, hjust = 0.5,
point_size = 0.1, alpha = 0.01, label=NULL, label_size = 3.95)+
scale_fill_brewer(palette = "Paired")+
scale_color_brewer(palette = "Paired")+
geom_point(aes(size = vsize, fill = Lineage_rel), alpha = 0.5,
colour="black", shape=21)+
guides(size = FALSE,
fill = guide_legend(title = "Phylum", override.aes = list(size = 5),
nrow = 4),
color = FALSE)+
theme(legend.position="bottom", legend.text=element_text(size=12),
text = element_text(size = 12),
plot.margin = unit(c(1,1,1,1), "cm"))+
scale_size(range = c(5, 15))+
geom_label_repel(aes(label = limno_labs),
fontface = 'bold', color = 'black',
box.padding = 0.35, point.padding = 0.5,
segment.color = 'black',
size = 4,
# Width of the line segments.
segment.size = 1.5,
# Draw an arrow from the label to the data point.
arrow = arrow(length = unit(0.015, 'npc')),
nudge_x = -0.1,
nudge_y = 0.6
)
print(p_16S_network)
2. Investigate MAG- and 16S-based abundances
2.1. From coverage info (anvi’o)
Formula used to calculate relative abundances: \[Relative\ abundance =100*(\frac{mean\ coverage * bin\ size}{read\ length*total\ sample\ reads })\]
# Normalize for bin sizes
data_total_abund <- data_total %>% group_by(bins) %>%
mutate(norm_mean_rel_abundance = mean_rel_abundance/(bin_size/1e6),
norm_upper_rel_abundance = upper_rel_abundance/(bin_size/1e6),
norm_lower_rel_abundance = lower_rel_abundance/(bin_size/1e6)
)
# Import bin name file used in manuscript
new_bin_names <- read.table("./anvio_output/rebin/general_bins_summary_selected_final.tsv", header = TRUE)[, c(2,3)]
data_total_abund <- left_join(data_total_abund,
new_bin_names, by = c("bins" = "bins"))
# Remove non-Limno bin
data_total_abund <- data_total_abund %>%
dplyr::filter(new_bin_name != "MAG.noLIM")
# Add putative limnohabitans lineage ID
data_total_abund <- dplyr::left_join(data_total_abund, MGT_df,
by = c("new_bin_name" = "Genome_ID"))
# Order names for phylogenomic tree
data_total_abund$new_bin_name <- as.character(data_total_abund$new_bin_name)
data_total_abund$new_bin_name <- factor(data_total_abund$new_bin_name,
levels = c(
"MAG5.SP-M110-DD",
"MAG2.FA-MLB-SN",
"MAG3.FA-MLB-SN",
"MAG4.FA-M110-DN",
"MAG1.FA-MLB-DN",
"MAG10.SU-M15-SN",
"MAG6.SP-M15-SD",
"MAG8.SU-M110-DCMD",
"MAG7.SU-MLB-SD",
"MAG9.SU-M15-SN"
)
)
# Plot abundance distributions of all bins
p_season1 <- ggplot(data = data_total_abund,
aes(x = new_bin_name, y = norm_mean_rel_abundance,
fill = Lineage))+
geom_point(size = 4, shape = 21, alpha = 0.7)+
geom_boxplot(alpha = 0.7, outlier.shape = NA)+
theme_bw()+
scale_fill_manual(values = c(rgb(red=t(col2rgb("#deebf7ff")),
maxColorValue = 255),
rgb(red=t(col2rgb("#c6dbefff")),
maxColorValue = 255),
rgb(red=t(col2rgb("#9ecae1ff")),
maxColorValue = 255),
rgb(red=t(col2rgb("#6baed6ff")),
maxColorValue = 255)
))+
# ylim(0,1)+
theme(axis.text=element_text(size=14), axis.title=element_text(size=20),
title=element_text(size=20), legend.text=element_text(size=14),
legend.background = element_rect(fill="transparent"),
axis.text.x = element_text(angle = 90, hjust = 1),
strip.text=element_text(size=18),
axis.line = element_line(size = 1, colour = "grey80"),
panel.border = element_blank())+
ylab("Normalized relative abundance (%)")+
scale_y_sqrt()+
xlab("")
p_season1
2.1. From mapped reads info (BBmap - pileup.sh)
From: https://doi.org/10.1186/s13059-015-0834-7
“The relative abundance of each MAG was estimated using the fraction of reads in each sample mapping to the respective MAG. Normalized on the size of that bin, this yielded the measure fraction of reads per nucleotide in bin. This measure was chosen since it is comparable across samples with varying sequencing output and different bin sizes. Using the CONCOCT input table, multiplying the average coverage per nucleotide with the length of the contig in question and summing over all contigs within a bin and within a sample gave the number of reads per bin within a sample. The fraction of reads in each sample mapping to each bin was then calculated by dividing this value with the total number of reads from each sample, after having removed duplicated reads.”
# Generated from SAM file with pileup.sh from BBmap toolset
df_map_reads <- read.delim("./anvio_output/collection_vizbin_pileup_cov.tsv",
stringsAsFactors = FALSE)
df_map_reads$Sample <- gsub(df_map_reads$Sample, pattern = "_cov.tsv",
replacement = "")
# This file was generated in anvio with anvi-export-collections
lab_map_reads <- read.delim("./anvio_output/collection_vizbin.tsv",
header = FALSE)
colnames(lab_map_reads) <- c("contig_split", "Bin")
# Select Limnohabitans bins
lab_map_reads <- lab_map_reads[lab_map_reads$Bin %in% data_total$bins, ]
lab_map_reads$contig <- paste(
do.call(rbind, strsplit(as.character(lab_map_reads$contig_split), "_"))[, 1],
do.call(rbind, strsplit(as.character(lab_map_reads$contig_split), "_"))[, 2],
sep = "_"
)
lab_map_reads <- lab_map_reads[ ,2:3] %>% distinct()
# Merge dataframes
df_map_merged <- left_join(lab_map_reads, df_map_reads, by = c("contig" = "X.ID")
)
df_map_merged[, 4:13] <- apply(df_map_merged[, 4:13], 2, FUN = function(x) as.numeric(x))
df_map_merged <- df_map_merged %>% group_by(Bin, Sample) %>% summarise(sum_map_read = sum(Plus_reads))
# Merge with metadata/bin sizes
meta$Sample <- gsub(".C", "", meta$Sample, fixed = TRUE)
meta$Sample <- gsub(".A", "", meta$Sample, fixed = TRUE)
df_map_merged$Sample <- gsub(".C", "", df_map_merged$Sample, fixed = TRUE)
df_map_merged <- left_join(df_map_merged, bin_size, by = c("Bin" = "bins"))
df_map_merged <- left_join(df_map_merged, meta, "Sample")
total_reads$sample <- gsub("_", ".", fixed = TRUE,total_reads$sample)
df_map_merged <- left_join(df_map_merged, total_reads, c("Sample" = "sample"))
df_map_merged <- df_map_merged %>% group_by(Bin, Sample) %>%
mutate(rel_abundance = 100*sum_map_read/Total_reads)
df_map_merged <- df_map_merged %>% group_by(Bin, Sample) %>%
mutate(rel_norm_abundance = 100*sum_map_read/Total_reads/(bin_size/1e6))
# Rename bin names to more sensible names
# Add extra column with new bin names
new_bin_names <- read.table("./anvio_output/rebin/general_bins_summary_selected_final.tsv", header = TRUE)[, c(2,3)]
df_map_merged <- left_join(df_map_merged, new_bin_names, by = c("Bin" = "bins"))
df_map_merged$new_bin_name <- as.character(df_map_merged$new_bin_name)
df_map_merged$new_bin_name <- factor(df_map_merged$new_bin_name, levels =
c("MAG1.FA-MLB-DN","MAG2.FA-MLB-SN",
"MAG3.FA-MLB-SN", "MAG4.FA-M110-DN",
"MAG5.SP-M110-DD","MAG6.SP-M15-SD",
"MAG7.SU-MLB-SD","MAG8.SU-M110-DCMD",
"MAG9.SU-M15-SN","MAG10.SU-M15-SN"))
df_map_merged$Site <- as.character(df_map_merged$Site)
df_map_merged$Site <- gsub("110", "Lake Michigan\nsite M110", df_map_merged$Site)
df_map_merged$Site <- gsub("15", "Lake Michigan\nsite M15", df_map_merged$Site)
df_map_merged$Site <- gsub("Buoy", "Muskegon Lake", df_map_merged$Site)
df_map_merged$Site <- factor(df_map_merged$Site, levels = c("Muskegon Lake",
"Lake Michigan\nsite M15",
"Lake Michigan\nsite M110"))
df_map_merged$Season <- as.character(df_map_merged$Season)
df_map_merged$Season <- factor(df_map_merged$Season, levels = c("Spring", "Summer","Fall"))
# Remove non-limno bin
df_map_merged <- df_map_merged %>% dplyr::filter(new_bin_name != "MAG.noLIM")
# Make plots
p_abs2 <- ggplot(data = df_map_merged, aes(x = new_bin_name, y = rel_norm_abundance, fill = new_bin_name))+
theme_bw()+
scale_fill_manual("", values = fill_palette)+
geom_jitter(size = 4, shape = 21, color = "black", alpha = 0.7, width = 0.15)+
theme(axis.text=element_text(size=14), axis.title=element_text(size=20),
title=element_text(size=20), legend.text=element_text(size=12),
legend.background = element_rect(fill="transparent"),
axis.text.x = element_blank(),
strip.text=element_text(size=14), legend.position = "bottom",
strip.background = element_rect(fill = adjustcolor("gray", 0.15)))+
ylab(paste0("Norm. relative abundance (%)"))+
xlab("")+
guides(fill=guide_legend(nrow = 3))+
facet_grid(Season~Site, scales ="free")+
scale_y_continuous(labels=scaleFUN, limits = c(0,1.25))+
coord_trans(y = "sqrt")
p_abs2
3. Expression analysis
One way: 1. Get total nr. of reads through the samtools flagstat command. 2. Get gene length through the DESMAN command: python $DESMAN/scripts/Lengths.py -i CDS.fa > CDS.len
Chosen way: 1. Map reads with bwa 2. Extract mapping statistics with pileup.sh from BBtools
Import and process data
The file needed to map genome_ids to gene_ids was generated with this code individually for each IMG genes.fna file and the concatenated with cat:
grep ">" 2757320404.genes.fna | sed "s/>//g" | awk '{ print $1, $2 }' > 2757320404.genes.ids
awk -F ' ' -v samname="2757320404" '{$(1)=samname FS $1;}1' OFS=" " 2757320404.genes.ids > 2757320404.genes.namesexpr_cov <- read.delim("./metaT/metaT_pileup.tsv", header = TRUE,
stringsAsFactors = FALSE)
# Merge gene annotations from all genomes in one file
file_list <- list.files(path = "./IMG_annotation", recursive = FALSE,
pattern = "IMG", full.names = TRUE)
# Import annotations of all MAG genes
merged_file <- merge_annotations(file_list[1:10], genoid_seqid = FALSE)## --- I will merge the annotation files from the following genomes:
## Genomes
## 1 ./IMG_annotation/IMG_2757320395
## 2 ./IMG_annotation/IMG_2757320396
## 3 ./IMG_annotation/IMG_2757320397
## 4 ./IMG_annotation/IMG_2757320398
## 5 ./IMG_annotation/IMG_2757320399
## 6 ./IMG_annotation/IMG_2757320400
## 7 ./IMG_annotation/IMG_2757320401
## 8 ./IMG_annotation/IMG_2757320402
## 9 ./IMG_annotation/IMG_2757320403
## 10 ./IMG_annotation/IMG_2757320404
## [1] 2322
## [1] 3720
## [1] 5259
## [1] 7208
## [1] 9745
## [1] 12405
## [1] 14858
## [1] 17644
## [1] 19813
## [1] 21848
## Wed Jun 06 20:51:29 2018 --- Sucessfully merged filesmerged_file$ko_id <- gsub(merged_file$ko_id, pattern = "KO:",
replacement = "")
# import KO annotation hierarchy
ko_df <- format_ko(path = "./Mapping_files/ko00000.keg")
# Annotated merged_file
merged_file_annot <- dplyr::left_join(merged_file, ko_df, by = "ko_id")
# Genes that were not annotated can be assigned a genome ID using another file
# made from the genes.fna IMG generated files
genome_gene_map <- read.table("IMG_annotation/genome_geneids.txt",
header = FALSE)[, 1:2]
colnames(genome_gene_map) <- c("Genome_ID", "gene_oid")
genome_gene_map <- data.frame(apply(genome_gene_map, 2, as.character))
# Add genome_ID to all genes
expr_cov <- dplyr::left_join(expr_cov, genome_gene_map, by = "gene_oid")
# Annotate this expression table with Kegg Orthology
# expr_cov <- dplyr::left_join(expr_cov, merged_file[, c(1,9)], by = c("gene_oid"))
expr_cov$Plus_reads <- as.integer(expr_cov$Plus_reads)
expr_cov$Length <- as.integer(expr_cov$Length)
expr_cov <- expr_cov[!is.na(expr_cov$Plus_reads),]
expr_cov_long <- dplyr::mutate(expr_cov,
reads_per_kb = Plus_reads/Length/1000
)
expr_cov_long <- expr_cov_long %>% group_by(Sample) %>%
mutate(RPK_scaling = sum(reads_per_kb)/1e6
)
expr_cov_long <- expr_cov_long %>%
mutate(TPM = reads_per_kb/RPK_scaling
)
# Now add the metadata to this long dataframe
# Metadata file
meta_metaT <- distinct(meta[, 2:nrow(meta)])
meta_metaT$Sample_ID <- gsub(meta_metaT$Sample_ID, pattern="_", replacement = ".")
rownames(meta_metaT) <- meta_metaT$Sample_ID
expr_cov_long <- left_join(expr_cov_long, meta_metaT[, 1:11],
by = c("Sample" = "Sample_ID"))
expr_cov_long$Genome_ID <- as.factor(expr_cov_long$Genome_ID)
# Remove duplicate rows
expr_cov_long <- expr_cov_long %>% distinct()Calculate dissimarility between expression profiles
colnames(expr_cov_long)[colnames(expr_cov_long) == "Plus_reads"] <- "mapped_reads"
# For each MAG calculat the dissimarilties in expression profiles between samples
for(i in 1:nlevels(expr_cov_long$Genome_ID)){
tmp_table <- expr_cov_long %>%
dplyr::filter(Genome_ID == levels(expr_cov_long$Genome_ID)[i]) %>%
dplyr::select(gene_oid, Sample, mapped_reads) %>%
tidyr::spread(Sample, mapped_reads)
r.geneoid <- tmp_table$gene_oid
min_sample_size <- min(colSums(tmp_table[,-1]))
tmp_table <- vegan::rrarefy(t(tmp_table[, -1]), min_sample_size)
tmp_table <- scale(tmp_table)
tmp_dist <- dist(tmp_table, method = "euclidean")
if(i == 1) results_diss <- data.frame(dist_m = matrix(tmp_dist),
Genome = levels(expr_cov_long$Genome_ID)[i]) else {
results_diss <- rbind(results_diss,
data.frame(dist_m = matrix(tmp_dist),
Genome = levels(expr_cov_long$Genome_ID)[i]))
}
remove(tmp_table)
}
# Merge with correct sample names
new_bin_names2 <- read.table("./anvio_output/rebin/general_bins_summary_selected_final.tsv", header = TRUE)[, c(3,8:10)];
new_bin_names2$IMG_taxID <- as.character(new_bin_names2$IMG_taxID)
results_diss <- left_join(results_diss, new_bin_names2,
by = c("Genome" = "IMG_taxID"))
results_diss$new_bin_name <- as.character(results_diss$new_bin_name)
results_diss$new_bin_name <- factor(results_diss$new_bin_name, levels =
c("MAG5.SP-M110-DD","MAG2.FA-MLB-SN",
"MAG3.FA-MLB-SN", "MAG4.FA-M110-DN",
"MAG1.FA-MLB-DN", "MAG10.SU-M15-SN",
"MAG6.SP-M15-SD", "MAG8.SU-M110-DCMD",
"MAG7.SU-MLB-SD","MAG9.SU-M15-SN"))
# Plot
results_diss %>%
ggplot(aes(x = new_bin_name, y = dist_m, fill = new_bin_name))+
# geom_jitter(size = 2.5,
# color = "black", shape = 21, width = 0.1, alpha = 0.5)+
geom_violin(alpha = 0.4, adjust = 1, draw_quantiles = TRUE)+
stat_summary(fun.data=mean_sdl, fun.args = list(mult = 1),
geom="pointrange", color="#333333", size = 1.05)+
# ylim(0,120)+
scale_fill_manual("",
values = c(rep(rgb(red=t(col2rgb("#deebf7ff")),
maxColorValue = 255), 2),
rep(rgb(red=t(col2rgb("#c6dbefff")),
maxColorValue = 255),5),
rgb(red=t(col2rgb("#9ecae1ff")), maxColorValue = 255),
rep(rgb(red=t(col2rgb("#6baed6ff")),
maxColorValue = 255),2)
))+
theme_bw()+
theme(axis.text.y=element_text(size=14), axis.title=element_text(size=20),
axis.text.x=element_text(size=14, angle = 45, hjust = 1),
title=element_text(size=16), legend.text=element_text(size=14),
legend.background = element_rect(fill="transparent"),
strip.text.x=element_text(size=18))+
guides(fill = FALSE)+
labs(x="", y = "Euclidean distance")
Run DESeq
Don’t forget:
In order to benefit from the default settings of the package, you should put the variable of interest at the end of the formula and make sure the control level is the first level.
Manual can be found here
# Format data for DESeq2
## Put count matrices in list including a count matrix for each bin
expr_cov_bins <- list()
for(i in 1:nlevels(expr_cov_long$Genome_ID)){
expr_cov_bins[[i]] <- expr_cov_long %>%
dplyr::filter(Genome_ID == levels(expr_cov_long$Genome_ID)[i]) %>%
dplyr::select(gene_oid, Sample, mapped_reads) %>%
tidyr::spread(Sample, mapped_reads)
r.bin <- expr_cov_bins[[i]]$gene_oid
expr_cov_bins[[i]] <- as.matrix(expr_cov_bins[[i]][, -1])
rownames(expr_cov_bins[[i]]) <- r.bin
}
# Check order of colnames in count and rownames in metadata matrix
# and make sure these are in the same order!
meta_metaT <- as.matrix(meta_metaT)
meta_metaT <- meta_metaT[match(colnames(expr_cov_bins[[1]]), rownames(meta_metaT)), ]
all(rownames(meta_metaT) %in% colnames(expr_cov_bins[[1]]))## [1] TRUE# Perform DESeq2 for differential abundance testing for each genome separately
sel_env <- meta_metaT[, "Site"] == "110" & meta_metaT[, "Depth"] != "Mid" & meta_metaT[, "Season"] != "Spring" & meta_metaT[, "Time"] == "Night"
expr_cov_MAG8[[1]] <- expr_cov_MAG8[[1]][, sel_env]## Error in eval(expr, envir, enclos): object 'expr_cov_MAG8' not foundmeta_metaT_subs <- meta_metaT[sel_env, ]
for(i in 1:nlevels(expr_cov_long$Genome_ID)){
}
## Season effect
General_deseq_results_season <- list()
deseq_comparisons_season <- list()
for(i in 1:nlevels(expr_cov_long$Genome_ID)){
cat(" --- Running DESeq2 on Genome_ID:",levels(expr_cov_long$Genome_ID)[i], "\n",sep = " ")
# Test for season but controlling for site
dds <- DESeq2::DESeqDataSetFromMatrix(countData = expr_cov_bins[[i]],
colData = meta_metaT,
design= ~ Site + Season)
# Run DESeq
dds <- DESeq2::DESeq(dds, quiet = TRUE)
# Calculate contrasts for all comparisons
comp1 <- DESeq2::results(dds, contrast=c("Season","Fall","Summer"))[order(results(dds)$padj), ]
comp2 <- DESeq2::results(dds, contrast=c("Season","Fall","Spring"))[order(results(dds)$padj), ]
comp3 <- DESeq2::results(dds, contrast=c("Season","Spring","Summer"))[order(results(dds)$padj), ]
# Store data in single list with dataframes
General_deseq_results_season[[i]] <- rbind(comp1, comp2, comp3)
deseq_comparisons_season[[i]] <- data.frame(comparison =
c(rep("Fall-Summer", nrow(comp1)),
rep("Fall-Spring", nrow(comp2)),
rep("Spring-Summer", nrow(comp3))),
design = "~ Site + Season"
)
}## --- Running DESeq2 on Genome_ID: 2757320395
## --- Running DESeq2 on Genome_ID: 2757320396
## --- Running DESeq2 on Genome_ID: 2757320397
## --- Running DESeq2 on Genome_ID: 2757320398
## --- Running DESeq2 on Genome_ID: 2757320399
## --- Running DESeq2 on Genome_ID: 2757320400
## --- Running DESeq2 on Genome_ID: 2757320401
## --- Running DESeq2 on Genome_ID: 2757320402
## --- Running DESeq2 on Genome_ID: 2757320403
## --- Running DESeq2 on Genome_ID: 2757320404## Site effect
General_deseq_results_site <- list()
deseq_comparisons_site <- list()
for(i in 1:nlevels(expr_cov_long$Genome_ID)){
cat(" --- Running DESeq2 on Genome_ID:",levels(expr_cov_long$Genome_ID)[i], "\n",sep = " ")
dds <- DESeq2::DESeqDataSetFromMatrix(countData = expr_cov_bins[[i]],
colData = meta_metaT,
design= ~ Season + Site) # Test for site but controlling for season
dds <- DESeq2::DESeq(dds, quiet = TRUE)
# Calculate contrasts for all comparisons
comp1 <- DESeq2::results(dds,
contrast=c("Site","Buoy","110"))[order(results(dds)$padj), ]
comp2 <- DESeq2::results(dds,
contrast=c("Site","Buoy","15"))[order(results(dds)$padj), ]
comp3 <- DESeq2::results(dds,
contrast=c("Site","15","110"))[order(results(dds)$padj), ]
# Store data in single list with dataframes
General_deseq_results_site[[i]] <- rbind(comp1, comp2, comp3)
deseq_comparisons_site[[i]] <- data.frame(comparison =
c(rep("Muskegon Lake\nVs\nM110",
nrow(comp1)),
rep("Muskegon Lake\nVs\nM15",
nrow(comp2)),
rep("M15\nVs\nM110",
nrow(comp3))),
design = "~ Season + Site"
)
}## --- Running DESeq2 on Genome_ID: 2757320395
## --- Running DESeq2 on Genome_ID: 2757320396
## --- Running DESeq2 on Genome_ID: 2757320397
## --- Running DESeq2 on Genome_ID: 2757320398
## --- Running DESeq2 on Genome_ID: 2757320399
## --- Running DESeq2 on Genome_ID: 2757320400
## --- Running DESeq2 on Genome_ID: 2757320401
## --- Running DESeq2 on Genome_ID: 2757320402
## --- Running DESeq2 on Genome_ID: 2757320403
## --- Running DESeq2 on Genome_ID: 2757320404## Depth effect
General_deseq_results_depth <- list()
deseq_comparisons_depth <- list()
for(i in 1:nlevels(expr_cov_long$Genome_ID)){
cat(" --- Running DESeq2 on Genome_ID:",levels(expr_cov_long$Genome_ID)[i], "\n",sep = " ")
# Test for Depth but controlling for site
dds <- DESeq2::DESeqDataSetFromMatrix(countData = expr_cov_bins[[i]],
colData = meta_metaT,
design= ~ Site + Depth)
# Run DESeq
dds <- DESeq2::DESeq(dds, quiet = TRUE)
# Calculate contrasts for all comparisons
comp1 <- DESeq2::results(dds, contrast=c("Depth","Surface","Deep"))[order(results(dds)$padj), ]
comp2 <- DESeq2::results(dds, contrast=c("Depth","Mid","Deep"))[order(results(dds)$padj), ]
comp3 <- DESeq2::results(dds, contrast=c("Depth","Surface","Mid"))[order(results(dds)$padj), ]
# Store data in single list with dataframes
General_deseq_results_depth[[i]] <- rbind(comp1, comp2, comp3)
deseq_comparisons_depth[[i]] <- data.frame(comparison =
c(rep("Surface-Deep", nrow(comp1)),
rep("Mid-Deep", nrow(comp2)),
rep("Surface-Mid", nrow(comp3))),
design = "~ Site + Depth"
)
}## --- Running DESeq2 on Genome_ID: 2757320395
## --- Running DESeq2 on Genome_ID: 2757320396
## --- Running DESeq2 on Genome_ID: 2757320397
## --- Running DESeq2 on Genome_ID: 2757320398
## --- Running DESeq2 on Genome_ID: 2757320399
## --- Running DESeq2 on Genome_ID: 2757320400
## --- Running DESeq2 on Genome_ID: 2757320401
## --- Running DESeq2 on Genome_ID: 2757320402
## --- Running DESeq2 on Genome_ID: 2757320403
## --- Running DESeq2 on Genome_ID: 2757320404# Pool results for each genome into a single result dataframe
res_deseq <- data.frame()
for(i in 1:length(General_deseq_results_season)){
tmp <- data.frame(gene_oid = General_deseq_results_season[[i]]@rownames,
baseMean = General_deseq_results_season[[i]]@listData$baseMean,
log2FoldChange = General_deseq_results_season[[i]]@listData$log2FoldChange,
pvalue = General_deseq_results_season[[i]]@listData$pvalue,
padj = General_deseq_results_season[[i]]@listData$padj,
Genome_ID = levels(expr_cov_long$Genome_ID)[i],
Comparison = deseq_comparisons_season[[i]]$comparison,
Design = deseq_comparisons_season[[i]]$design)
if(i == 1) res_deseq <- tmp else{
res_deseq <- rbind(res_deseq, tmp)
}
}
for(i in 1:length(General_deseq_results_site)){
tmp <- data.frame(gene_oid = General_deseq_results_site[[i]]@rownames,
baseMean = General_deseq_results_site[[i]]@listData$baseMean,
log2FoldChange = General_deseq_results_site[[i]]@listData$log2FoldChange,
pvalue = General_deseq_results_site[[i]]@listData$pvalue,
padj = General_deseq_results_site[[i]]@listData$padj,
Genome_ID = levels(expr_cov_long$Genome_ID)[i],
Comparison = deseq_comparisons_site[[i]]$comparison,
Design = deseq_comparisons_site[[i]]$design)
res_deseq <- rbind(res_deseq, tmp)
}
for(i in 1:length(General_deseq_results_depth)){
tmp <- data.frame(gene_oid = General_deseq_results_depth[[i]]@rownames,
baseMean = General_deseq_results_depth[[i]]@listData$baseMean,
log2FoldChange = General_deseq_results_depth[[i]]@listData$log2FoldChange,
pvalue = General_deseq_results_depth[[i]]@listData$pvalue,
padj = General_deseq_results_depth[[i]]@listData$padj,
Genome_ID = levels(expr_cov_long$Genome_ID)[i],
Comparison = deseq_comparisons_depth[[i]]$comparison,
Design = deseq_comparisons_depth[[i]]$design)
res_deseq <- rbind(res_deseq, tmp)
}
# Only retain significantly differentially expressed genes
res_deseq <- res_deseq %>% dplyr::filter(padj < 0.01)
new_bin_names2 <- read.table("./anvio_output/rebin/general_bins_summary_selected_final.tsv", header = TRUE)[, c(3,8:10)]; new_bin_names2$IMG_taxID <- as.character(new_bin_names2$IMG_taxID)
res_deseq <- left_join(res_deseq, new_bin_names2, by = c("Genome_ID" = "IMG_taxID"))
res_deseq$new_bin_name <- as.character(res_deseq$new_bin_name)
res_deseq$new_bin_name <- factor(res_deseq$new_bin_name, levels =
c("MAG1.FA-MLB-DN","MAG2.FA-MLB-SN",
"MAG3.FA-MLB-SN", "MAG4.FA-M110-DN",
"MAG5.SP-M110-DD", "MAG6.SP-M15-SD",
"MAG7.SU-MLB-SD", "MAG8.SU-M110-DCMD",
"MAG9.SU-M15-SN", "MAG10.SU-M15-SN"))
res_deseq$regulation <- res_deseq$log2FoldChange > 0
res_deseq$regulation[res_deseq$regulation == "TRUE"] <- "upregulation"
res_deseq$regulation[res_deseq$regulation == "FALSE"] <- "downregulation"
# Add KO/COG/pfam/tigrfam gene annotations to DeSEQ2 results
res_deseq_anott <- dplyr::left_join(res_deseq, merged_file_annot, by = c("gene_oid"))
metaT_reads_mapped <- expr_cov_long %>%
group_by(Sample, Genome_ID) %>%
dplyr::summarise(sum_mapped_reads = sum(mapped_reads))
metaT_reads_mapped <- left_join(metaT_reads_mapped, new_bin_names2, by = c("Genome_ID" = "IMG_taxID"))
metaT_reads_mapped$new_bin_name <- as.character(metaT_reads_mapped$new_bin_name)
metaT_reads_mapped$new_bin_name <- factor(metaT_reads_mapped$new_bin_name,
levels =
c("MAG1.FA-MLB-DN","MAG2.FA-MLB-SN",
"MAG3.FA-MLB-SN", "MAG4.FA-M110-DN",
"MAG5.SP-M110-DD", "MAG6.SP-M15-SD",
"MAG7.SU-MLB-SD", "MAG8.SU-M110-DCMD",
"MAG9.SU-M15-SN", "MAG10.SU-M15-SN"))
p_deseq_reads <- ggplot2::ggplot(metaT_reads_mapped, aes(x = new_bin_name, y = sum_mapped_reads, fill = new_bin_name))+
geom_bar(color = "black", stat = "identity")+
theme_bw()+
scale_fill_brewer(palette = "Paired")+
theme(axis.text.x = element_blank(),
axis.text.y = element_text(size = 14),
legend.title = element_blank(),
axis.title = element_text(size = 14),
strip.text = element_text(size = 12),
axis.title.x = element_blank())+
ylab("Number of mapped reads")+
xlab("")+
facet_wrap(~Sample, nrow = 4)
print(p_deseq_reads)
selected_ko <- c("Carbohydrate metabolism", "Energy metabolism",
"Amino acid metabolism", "Membrane transport",
"Replication and repair", "Translation",
"Energy metabolism")
# for(genome in levels(unique(res_deseq_anott$new_bin_name))){
p_deseq_1 <- res_deseq_anott %>% dplyr::filter(!is.na(ko_id) & ko_level_B %in% selected_ko & Design == "~ Season + Site") %>%
ggplot2::ggplot(aes(x = ko_level_C, fill = ko_level_B))+
geom_bar(color = "black")+
theme_bw()+
scale_fill_brewer("", palette = "Paired")+
theme(axis.text.x = element_text(size = 10.5, angle = 45, hjust =1),
axis.text.y = element_text(size = 14),
legend.title = element_blank(),
axis.title = element_text(size = 14),
strip.text = element_text(size = 14)) +
ylab("Number of genes")+
xlab("")+
facet_wrap(~new_bin_name, ncol = 2)+
guides(fill = FALSE)+
ggtitle("Muskegon lake vs. M110")
print(p_deseq_1)
# }
# Plot for each genome the number of differentially expressed genes according to:
# season and site
# p_deseq_1 <- ggplot2::ggplot(res_deseq, aes(x = new_bin_name, fill = new_bin_name))+
# geom_bar(color = "black")+
# theme_bw()+
# scale_fill_brewer(palette = "Paired")+
# theme(axis.text.x = element_blank(),
# axis.text.y = element_text(size = 14),
# legend.title = element_blank(),
# axis.title = element_text(size = 14),
# strip.text = element_text(size = 14)) +
# ylab("Number of mapped reads")+
# xlab("")+
# facet_grid(Design~Comparison)
#
# print(p_deseq_1)
# p_deseq_2 <- ggplot2::ggplot(res_deseq, aes(x = new_bin_name, shape = regulation,
# y = log2FoldChange))+
# # geom_jitter(size = 3, width = 0.25, shape = 21)+
# geom_boxplot(alpha= 1, aes(fill = new_bin_name))+
# theme_bw()+
# scale_fill_brewer(palette = "Paired")+
# theme(axis.text.x = element_blank(),
# axis.text.y = element_text(size = 14),
# legend.title = element_blank(),
# axis.title = element_text(size = 14),
# strip.text = element_text(size = 14))+
# ylab("log2FoldChange")+
# xlab("")+
# facet_grid(Design~Comparison)+
# guides(shape = FALSE)
#
# print(p_deseq_2)Run MAG8-DESeq
# Select MAG8 genome and temperature gradient
expr_cov_MAG8 <- expr_cov_bins[levels(expr_cov_long$Genome_ID) == "2757320398"]
sel_MAG8 <- meta_metaT[, "Site"] == "110" & meta_metaT[, "Depth"] != "Mid" & meta_metaT[, "Season"] != "Spring" & meta_metaT[, "Time"] == "Night"
expr_cov_MAG8[[1]] <- expr_cov_MAG8[[1]][, sel_MAG8]
meta_metaT_MAG8 <- meta_metaT[sel_MAG8, ]
## Depth (temperature effect) effect and control for Seasonal variation
dds <- DESeq2::DESeqDataSetFromMatrix(countData = expr_cov_MAG8[[1]],
colData = meta_metaT_MAG8,
design= ~ Season + Depth)
# Run DESeq
dds <- DESeq2::DESeq(dds, quiet = TRUE)
# Calculate contrasts for all comparisons (deep vs surface logFC)
comp1 <- DESeq2::results(dds, contrast=c("Depth", "Deep", "Surface"))[order(results(dds)$padj), ]
# comp2 <- DESeq2::results(dds, contrast=c("Time.Depth", "Night.Deep", "Night.Surface"))[order(results(dds)$padj), ]
# Store data in single dataframes
MAG8_deseq_results_depth_comp1 <- data.frame(gene_oid = comp1@rownames,
baseMean = comp1@listData$baseMean,
log2FoldChange = comp1@listData$log2FoldChange,
pvalue = comp1@listData$pvalue,
padj = comp1@listData$padj,
Genome_ID = "2757320398",
Comparison = "Day.Deep - Day.Surface",
Time_of_day = "Day"
)
# MAG8_deseq_results_depth_comp2 <- data.frame(gene_oid = comp2@rownames,
# baseMean = comp2@listData$baseMean,
# log2FoldChange = comp2@listData$log2FoldChange,
# pvalue = comp2@listData$pvalue,
# padj = comp2@listData$padj,
# Genome_ID = "2757320398",
# Comparison = "Night.Deep - Night.Surface",
# Time_of_day = "Night"
# )
# MAG8_deseq_results_depth <- rbind(MAG8_deseq_results_depth_comp1,
# MAG8_deseq_results_depth_comp2)
MAG8_deseq_results_depth <- MAG8_deseq_results_depth_comp1
# Filter at p < 0.01
MAG8_deseq_results_depth <- MAG8_deseq_results_depth %>% dplyr::filter(padj < 0.01)# Select MAG8 genome and temperature gradient in spring M110
expr_cov_MAG8 <- expr_cov_bins[levels(expr_cov_long$Genome_ID) == "2757320398"]
sel_MAG8 <- meta_metaT[, "Site"] == "110" & meta_metaT[, "Depth"] == "Surface" & meta_metaT[, "Season"] %in% c("Spring", "Summer")
expr_cov_MAG8[[1]] <- expr_cov_MAG8[[1]][, sel_MAG8]
meta_metaT_MAG8 <- meta_metaT[sel_MAG8, ]
## Depth (temperature effect) effect and control for Seasonal variation
dds <- DESeq2::DESeqDataSetFromMatrix(countData = expr_cov_MAG8[[1]],
colData = meta_metaT_MAG8,
design= ~ Time + Season)
# Run DESeq
dds <- DESeq2::DESeq(dds, quiet = TRUE)
# Calculate contrasts for all comparisons (deep vs surface logFC)
comp1 <- DESeq2::results(dds, contrast=c("Season", "Spring", "Summer"))[order(results(dds)$padj), ]
# comp2 <- DESeq2::results(dds, contrast=c("Time.Depth", "Night.Deep", "Night.Surface"))[order(results(dds)$padj), ]
# Store data in single dataframes
MAG8_deseq_results_temp_comp1 <- data.frame(gene_oid = comp1@rownames,
baseMean = comp1@listData$baseMean,
log2FoldChange = comp1@listData$log2FoldChange,
pvalue = comp1@listData$pvalue,
padj = comp1@listData$padj,
Genome_ID = "2757320398",
Comparison = "Cold(Spring) - Warm(Summer)"
)
# MAG8_deseq_results_depth_comp2 <- data.frame(gene_oid = comp2@rownames,
# baseMean = comp2@listData$baseMean,
# log2FoldChange = comp2@listData$log2FoldChange,
# pvalue = comp2@listData$pvalue,
# padj = comp2@listData$padj,
# Genome_ID = "2757320398",
# Comparison = "Night.Deep - Night.Surface",
# Time_of_day = "Night"
# )
# MAG8_deseq_results_depth <- rbind(MAG8_deseq_results_depth_comp1,
# MAG8_deseq_results_depth_comp2)
MAG8_deseq_results_temp <- MAG8_deseq_results_temp_comp1
# Filter at p < 0.01
MAG8_deseq_results_temp <- MAG8_deseq_results_temp %>% dplyr::filter(padj < 0.01)# Merge with Kegg and COG annotations
# Import COG mapping file
cogid_2_cogcat <- read.csv("./Mapping_files/cogid_2_cogcat.csv", sep = ",", header = FALSE, fill = TRUE,col.names = c("COG_ID", "COG_class", "function"))[, 1:2]
cogid_2_cogcat <- cogid_2_cogcat[(cogid_2_cogcat$COG_class)!="", ]
cogid_2_cogcat <- droplevels(cogid_2_cogcat)
# Read COG category file
cog_categories <- read.table("./Mapping_files/cog_categories.tsv", header = TRUE, sep = "\t")
# Merge COG metadata
cog_meta <- dplyr::left_join(cog_categories, cogid_2_cogcat, by = c("COG_class" = "COG_class"))
cog_meta <- droplevels(cog_meta)
# Import KO hierarchical data
ko_path_df <- format_ko(path = "./Mapping_files/ko00000.keg")
# Import COG and KEGG annotation data of MAG8
MAG8_KO <- read.table("./IMG_annotation/IMG_2757320398/IMG Data/2757320398/2757320398.ko.tab.txt",
header = TRUE, fill = TRUE, sep = "\t" )
MAG8_KO$ko_id <- gsub("KO:", "", MAG8_KO$ko_id)
MAG8_COG <- read.table("./IMG_annotation/IMG_2757320398/IMG Data/2757320398/2757320398.cog.tab.txt",
stringsAsFactors = FALSE,
sep = "\t", header = TRUE,
quote=NULL, comment='')
# add hierarchy to cog and KO annotations
MAG8_KO <- left_join(MAG8_KO, ko_path_df, by = "ko_id")
MAG8_COG <- left_join(MAG8_COG, cog_meta, by = c("cog_id" = "COG_ID"))
# Merge annotations and remove duplicate columns
MAG8_annot <- full_join(MAG8_KO[, c(1,10:11, 14:18)],
MAG8_COG[, c(1,10:11, 13:15)], by = "gene_oid")
MAG8_annot$gene_oid <- as.character(MAG8_annot$gene_oid)
# Annotate differentially expressed genes
MAG8_deseq_results_depth_fin <- left_join(MAG8_deseq_results_depth, MAG8_annot,
by = c("gene_oid") ) %>% distinct()
MAG8_deseq_results_temp_fin <- left_join(MAG8_deseq_results_temp, MAG8_annot,
by = c("gene_oid") ) %>% distinct()
# Change NA to unknown
MAG8_deseq_results_depth_fin[,9:20] <- apply(MAG8_deseq_results_depth_fin[,9:20],
2, function(x) as.character(x))
MAG8_deseq_results_temp_fin[,8:19] <- apply(MAG8_deseq_results_temp_fin[,8:19],
2, function(x) as.character(x))
MAG8_deseq_results_depth_fin[is.na(MAG8_deseq_results_depth_fin)] <- "Unknown"
MAG8_deseq_results_temp_fin[is.na(MAG8_deseq_results_temp_fin)] <- "Unknown"
# Add label for up or downregulation
MAG8_deseq_results_depth_fin$regulation <- MAG8_deseq_results_depth_fin$log2FoldChange > 0
MAG8_deseq_results_depth_fin$regulation[MAG8_deseq_results_depth_fin$regulation == TRUE] <- "Upregulated"
MAG8_deseq_results_depth_fin$regulation[MAG8_deseq_results_depth_fin$regulation == FALSE] <- "Downregulated"
MAG8_deseq_results_depth_fin$regulation <- factor(MAG8_deseq_results_depth_fin$regulation,
levels =
c("Upregulated", "Downregulated"))
MAG8_deseq_results_temp_fin$regulation <- MAG8_deseq_results_temp_fin$log2FoldChange > 0
MAG8_deseq_results_temp_fin$regulation[MAG8_deseq_results_temp_fin$regulation == TRUE] <- "Upregulated"
MAG8_deseq_results_temp_fin$regulation[MAG8_deseq_results_temp_fin$regulation == FALSE] <- "Downregulated"
MAG8_deseq_results_temp_fin$regulation <- factor(MAG8_deseq_results_temp_fin$regulation,
levels =
c("Upregulated", "Downregulated"))# Overview
## Log fold changes deep vs surface
p_mag8_deseq_ov <- MAG8_deseq_results_depth_fin %>%
dplyr::select(gene_oid:Time_of_day) %>%
distinct() %>%
ggplot(aes(x = gene_oid, y = log2FoldChange, fill = log2FoldChange))+
geom_bar(size = 0.1, stat = "identity", color = "black")+
# geom_boxplot(alpha = 0.5)+
scale_fill_distiller(palette = "RdBu", name = "Fold Change",
limits = c(-15,15), oob=squish,
direction = -1) +
theme_bw()+
facet_grid(Time_of_day~.)+
theme(axis.text.x = element_blank())+
labs(x = "", y = "Fold Change (log2)")
print(p_mag8_deseq_ov)
# Focus on COG annotation and pool per level
p_mag8_deseq_cog <- MAG8_deseq_results_depth_fin %>% dplyr::select(gene_oid:Time_of_day, contains("cog")) %>%
distinct() %>%
ggplot(aes(x = COG_functional_category, y = log2FoldChange, fill = log2FoldChange))+
geom_point(size = 3, color = "black", shape = 21)+
# geom_boxplot(alpha = 0.5)+
scale_fill_distiller(palette = "RdBu", name = "Fold Change",
limits = c(-15,15), oob=squish,
direction = -1) +
theme_bw()+
facet_grid(Time_of_day~.)+
theme(axis.text.x = element_text(size = 12, angle = 90, hjust = 1))+
labs(x = "", y = "Fold Change (log2)")
print(p_mag8_deseq_cog)
p_mag8_deseq_temp_cog <- MAG8_deseq_results_temp_fin %>% dplyr::select(gene_oid:Comparison, contains("cog")) %>%
distinct() %>%
ggplot(aes(x = COG_functional_category, y = log2FoldChange, fill = log2FoldChange))+
geom_point(size = 3, color = "black", shape = 21)+
# geom_boxplot(alpha = 0.5)+
scale_fill_distiller(palette = "RdBu", name = "Fold Change",
limits = c(-15,15), oob=squish,
direction = -1) +
theme_bw()+
theme(axis.text.x = element_text(size = 12, angle = 90, hjust = 1))+
labs(x = "", y = "Fold Change (log2)")
print(p_mag8_deseq_temp_cog)
# Focus on KO annotation and pool per level
p_mag8_deseq_KO <- MAG8_deseq_results_depth_fin %>%
dplyr::select(gene_oid:Time_of_day, contains("ko")) %>%
distinct() %>%
ggplot(aes(x = ko_level_B, y = log2FoldChange, fill = log2FoldChange))+
geom_point(size = 3, color = "black", shape = 21)+
# geom_boxplot(alpha = 0.5)+
scale_fill_distiller(palette = "RdBu", name = "Fold Change",
limits = c(-15,15), oob=squish,
direction = -1) +
theme_bw()+
facet_grid(Time_of_day~.)+
theme(axis.text.x = element_text(size = 12, angle = 90, hjust = 1))+
labs(x = "", y = "Fold Change (log2)")
print(p_mag8_deseq_KO)
Plot
# Test for enrichment of functional categories in differentially expressed gene pool
## Test for enrichment of KO level C terms in transcriptome
bg_gsea <- MAG8_annot %>%
dplyr::select(ko_level_C, gene_oid) %>%
distinct()
bg_gsea[is.na(bg_gsea)] <- "Unknown"
metaT_gsea <- enricher(gene = unique(MAG8_deseq_results_depth_fin$gene_oid),
universe = bg_gsea$gene_oid,
TERM2GENE = bg_gsea,
pvalueCutoff = 0.05,
qvalueCutoff = 0.2)
metaT_gsea_KO <- data.frame(metaT_gsea@result)
metaT_gsea_KO %>%
dplyr::select(Description, GeneRatio, BgRatio, p.adjust, qvalue, Count)%>%
print()## Description
## Quorum sensing Quorum sensing
## Transporters Transporters
## ABC transporters ABC transporters
## Two-component system Two-component system
## Porphyrin and chlorophyll metabolism Porphyrin and chlorophyll metabolism
## GeneRatio BgRatio p.adjust
## Quorum sensing 14/108 36/1467 1.036930e-06
## Transporters 23/108 96/1467 1.036930e-06
## ABC transporters 16/108 50/1467 1.103524e-06
## Two-component system 6/108 20/1467 1.260791e-02
## Porphyrin and chlorophyll metabolism 4/108 10/1467 1.816558e-02
## qvalue Count
## Quorum sensing 8.434356e-07 14
## Transporters 8.434356e-07 23
## ABC transporters 8.976037e-07 16
## Two-component system 1.025524e-02 6
## Porphyrin and chlorophyll metabolism 1.477583e-02 4## Test for enrichment of COG functional categories
bg_gsea <- MAG8_annot %>%
dplyr::select(COG_functional_category, gene_oid) %>%
distinct()
bg_gsea <- apply(bg_gsea, 2, function(x) as.character(x))
bg_gsea[is.na(bg_gsea)] <- "Unknown"
bg_gsea <- data.frame(bg_gsea)
metaT_gsea <- enricher(gene = unique(MAG8_deseq_results_depth_fin$gene_oid),
universe = bg_gsea$gene_oid,
TERM2GENE = bg_gsea,
pvalueCutoff = 0.05,
qvalueCutoff = 0.2)
metaT_gsea_COG <- data.frame(metaT_gsea@result)
metaT_gsea_COG %>%
dplyr::select(Description, GeneRatio, BgRatio, p.adjust, qvalue, Count)%>%
print()## [1] Description GeneRatio BgRatio p.adjust qvalue Count
## <0 rows> (or 0-length row.names)# Focus on enriched KO level
p_mag8_deseq_depth_KO_gsea <- MAG8_deseq_results_depth_fin %>%
dplyr::select(regulation, gene_oid:Time_of_day, contains("ko")) %>%
distinct() %>%
dplyr::filter(ko_level_C %in% metaT_gsea_KO$Description) %>%
ggplot(aes(x = ko_level_C, y = log2FoldChange, fill = regulation))+
geom_point(shape = 21, size = 3, position=position_jitterdodge(dodge.width=0.9,
jitter.width=0.25)) +
geom_boxplot(outlier.colour = NA, width = 0.5,
position = position_dodge(width=0.9), alpha = 0.3)+
# geom_boxplot(alpha = 0.5)+
scale_fill_manual(values = brewer.pal(11, "RdBu")[c(2,10)], name = "Fold Change\n") +
theme_bw()+
theme(axis.text.x = element_text(size = 12, angle = 45, hjust = 1),
axis.text.y = element_text(size = 12),
strip.text = element_text(size = 14),
legend.text = element_text(size = 12),
legend.title = element_blank(),
legend.position = "top")+
labs(x = "", y = "Fold Change (log2)")+
ylim(-7, 7)
print(p_mag8_deseq_depth_KO_gsea)
# Focus on enriched KO level in sprin-summer samples too
p_mag8_deseq_temp_KO_gsea <- MAG8_deseq_results_temp_fin %>%
dplyr::select(regulation, gene_oid:Comparison, contains("ko")) %>%
distinct() %>%
dplyr::filter(ko_level_C %in% metaT_gsea_KO$Description) %>%
ggplot(aes(x = ko_level_C, y = log2FoldChange, fill = regulation))+
geom_point(shape = 21, size = 3, position=position_jitterdodge(dodge.width=0.9,
jitter.width=0.25)) +
geom_boxplot(outlier.colour = NA, width = 0.5,
position = position_dodge(width=0.9), alpha = 0.3)+
# geom_boxplot(alpha = 0.5)+
scale_fill_manual(values = brewer.pal(11, "RdBu")[c(2,10)], name = "Fold Change\n") +
theme_bw()+
theme(axis.text.x = element_text(size = 12, angle = 45, hjust = 1),
axis.text.y = element_text(size = 12),
strip.text = element_text(size = 14),
legend.text = element_text(size = 12),
legend.title = element_blank(),
legend.position = "top")+
labs(x = "", y = "Fold Change (log2)")+
ylim(-7, 7)
print(p_mag8_deseq_temp_KO_gsea)
# Spring - Summer samples at M110 - different temperatures but same depth
# Test for enrichment of functional categories in differentially expressed gene pool
## Test for enrichment of KO level C terms in transcriptome
bg_gsea <- MAG8_annot %>%
dplyr::select(ko_level_C, gene_oid) %>%
distinct()
bg_gsea[is.na(bg_gsea)] <- "Unknown"
metaT_gsea <- enricher(gene = unique(MAG8_deseq_results_temp_fin$gene_oid),
universe = bg_gsea$gene_oid,
TERM2GENE = bg_gsea,
pvalueCutoff = 0.05,
qvalueCutoff = 0.2)
metaT_gsea_KO <- data.frame(metaT_gsea@result)
metaT_gsea_KO %>%
dplyr::select(Description, GeneRatio, BgRatio, p.adjust, qvalue, Count)%>%
print()## [1] Description GeneRatio BgRatio p.adjust qvalue Count
## <0 rows> (or 0-length row.names)## Test for enrichment of COG functional categories
bg_gsea <- MAG8_annot %>%
dplyr::select(COG_functional_category, gene_oid) %>%
distinct()
bg_gsea <- apply(bg_gsea, 2, function(x) as.character(x))
bg_gsea[is.na(bg_gsea)] <- "Unknown"
bg_gsea <- data.frame(bg_gsea)
metaT_gsea <- enricher(gene = unique(MAG8_deseq_results_temp_fin$gene_oid),
universe = bg_gsea$gene_oid,
TERM2GENE = bg_gsea,
pvalueCutoff = 0.05,
qvalueCutoff = 0.2)
metaT_gsea_COG <- data.frame(metaT_gsea@result)
metaT_gsea_COG %>%
dplyr::select(Description, GeneRatio, BgRatio, p.adjust, qvalue, Count)%>%
print()## [1] Description GeneRatio BgRatio p.adjust qvalue Count
## <0 rows> (or 0-length row.names)# # Focus on enriched KO level
# p_mag8_deseq_KO_gsea <- MAG8_deseq_results_temp_fin %>%
# dplyr::select(regulation, gene_oid:Comparison, contains("ko")) %>%
# distinct() %>%
# dplyr::filter(ko_level_C %in% metaT_gsea_KO$Description) %>%
# ggplot(aes(x = ko_level_C, y = log2FoldChange, fill = regulation))+
# geom_point(shape = 21, size = 3, position=position_jitterdodge(dodge.width=0.9,
# jitter.width=0.25)) +
# geom_boxplot(outlier.colour = NA, width = 0.5,
# position = position_dodge(width=0.9), alpha = 0.3)+
# # geom_boxplot(alpha = 0.5)+
# scale_fill_manual(values = brewer.pal(11, "RdBu")[c(2,10)], name = "Fold Change\n") +
# theme_bw()+
# theme(axis.text.x = element_text(size = 12, angle = 45, hjust = 1),
# axis.text.y = element_text(size = 12),
# strip.text = element_text(size = 14),
# legend.text = element_text(size = 12),
# legend.title = element_blank(),
# legend.position = "top")+
# labs(x = "", y = "Fold Change (log2)")+
# ylim(-7, 7)
#
# print(p_mag8_deseq_KO_gsea)7. Sequence discrete populations
In order to check unspecific mapping sample reads were mapped consecutively to every bin using BBmap.sh. This approach allows us to check the mapping specificity by evaluating the distribution of read identity to the putative genome bin.
Competitive mapping was performed through blastn searches against a 1M and 10M read subsample from the interleaved fasta generated after QC (seqtk). Shell script used to achieve this:
#!/bin/bash
set -e
for file in `cat map.list`
do
echo $file
db=/scratch/vdenef_fluxm/rprops/DESMAN/metaG/vizbin_rebin_anvio_v230/SEQ_discrete/contigs/merged_bins-fixed.db
seqtk sample -s 777 /scratch/vdenef_fluxm/rprops/DESMAN/metaG/data/${file}/dt_int.fasta 1000000 > /scratch/vdenef_fluxm/rprops/DESMAN/metaG/data/${file}/dt_int_subs_1000000.fasta
qseqs=/scratch/vdenef_fluxm/rprops/DESMAN/metaG/data/${file}/dt_int_subs_1000000.fasta
blastn -query ${qseqs} -task megablast -db ${db} -out ./blast_output/${file}_blast.tsv -outfmt 6 -max_target_seqs 1 -num_threads 40 -perc_identity 60
rm ${qseqs}
doneCompetitive mapping using blast with 1M reads of interleaved fasta.
blast_df <- read.table("./SEQs_discrete/merged_blast_1M.tsv", header = FALSE)
colnames(blast_df) <- c("Sample", "Contig", "Identity")
blast_df_map <- read.table("./SEQs_discrete/merged_contig_list.tsv", header = FALSE)
colnames(blast_df_map) <- c("m_contig", "o_contig", "bin")
# Round identity to integer
# blast_df$Identity <- round(blast_df$Identity, 1)
# blast_df$Identity <- factor(blast_df$Identity)
blast_df <- left_join(blast_df, blast_df_map, by = c("Contig" = "m_contig"))
# Bin the %Identity in intervals of 0.5%
blast_df_sum <- transform(blast_df, bin_group=cut(Identity, breaks=seq(0, 100, 0.5)))
# Format sample names and convert to data.frame
blast_df_sum$Sample <- gsub("_blast.tsv", "", blast_df_sum$Sample)
blast_df_sum <- data.frame(blast_df_sum)
# blast_df_sum$Identity <- as.numeric(as.character(blast_df_sum$Identity))
blast_df_sum$bin <- gsub(".fa","",blast_df_sum$bin)
blast_df_sum <- dplyr::left_join(blast_df_sum, bin_size, by = c("bin" = "bins"))
# Add extra column that converts the binning range to a numeric x-coordinate that
# is positioned in the middle of the binning interval
blast_df_sum$bin_group <- gsub("\\(|]", "",
as.character(blast_df_sum$bin_group))
blast_df_sum$bin_xcoord <- as.numeric(do.call(rbind,
strsplit(blast_df_sum$bin_group,
","))[,1])+0.25
# Normalize mapped reads per sample based on sample reads
# blast_subs <- 1e6
# blast_df_sum <- blast_df_sum %>% ungroup %>% group_by(Sample) %>%
# mutate(n_prop = 100*n/blast_subs)
# Add season variable
blast_df_sum$season <- "Summer"
blast_df_sum$season[grep("Fa", blast_df_sum$Sample)] <- "Fall"
blast_df_sum$season[grep("Su", blast_df_sum$Sample)] <- "Summer"
blast_df_sum$season[grep("Sp", blast_df_sum$Sample)] <- "Spring"
# Reformat sample names
blast_df_sum$Sample <- gsub(".C", "", blast_df_sum$Sample, fixed = TRUE)
blast_df_sum$Sample <- gsub(".A", "", blast_df_sum$Sample, fixed = TRUE)
blast_df_sum$Sample <- gsub(".", "_", blast_df_sum$Sample, fixed = TRUE)
# Add metadata to dataframe
meta_blast <- meta[, -1] %>% distinct()
blast_df_sum <- dplyr::left_join(blast_df_sum, meta_blast, by = c("Sample" = "Sample_ID"))
# Reorder site factor
blast_df_sum$Site <- as.character(blast_df_sum$Site)
blast_df_sum$Site <- gsub("Buoy","Muskegon Lake", blast_df_sum$Site)
blast_df_sum$Site <- gsub("110","Lake Michigan\nsite M110", blast_df_sum$Site)
blast_df_sum$Site <- gsub("15","Lake Michigan\nsite M15", blast_df_sum$Site)
blast_df_sum$Site <- factor(blast_df_sum$Site, levels = c("Muskegon Lake",
"Lake Michigan\nsite M15",
"Lake Michigan\nsite M110"))
blast_df_sum$Depth <- as.character(blast_df_sum$Depth)
blast_df_sum$Depth <- factor(blast_df_sum$Depth, levels = c("Surface", "Mid", "Deep"))
blast_df_sum$season <- as.character(blast_df_sum$season)
blast_df_sum$season <- factor(blast_df_sum$season, levels = c("Spring", "Summer", "Fall"))
# remove non-Limnohabbitans bin
blast_df_sum <- blast_df_sum %>% dplyr::filter(bin != "B2_Fa13.BD.MLB.DN_rebin10")
# Add extra column with new bin names
new_bin_names <- read.table("./anvio_output/rebin/general_bins_summary_selected_final.tsv", header = TRUE)[, c(2,3)]
blast_df_sum <- left_join(blast_df_sum, new_bin_names, by = c("bin" = "bins"))
blast_df_sum$new_bin_name <- as.character(blast_df_sum$new_bin_name)
blast_df_sum$new_bin_name <- factor(blast_df_sum$new_bin_name, levels =
c("MAG1.FA-MLB-DN","MAG2.FA-MLB-SN",
"MAG3.FA-MLB-SN", "MAG4.FA-M110-DN",
"MAG5.SP-M110-DD","MAG6.SP-M15-SD",
"MAG7.SU-MLB-SD","MAG8.SU-M110-DCMD",
"MAG9.SU-M15-SN","MAG10.SU-M15-SN"))
# plot for individual bins
for(bin2plot in unique(blast_df_sum$new_bin_name)){
p_blast_sdisc <- blast_df_sum %>% dplyr::filter(new_bin_name == bin2plot) %>%
ggplot(aes(x = bin_xcoord, ..scaled.., fill = season))+
theme_bw()+
scale_fill_brewer(palette = "Accent")+
facet_wrap(~Sample, nrow = 4)+
geom_density(color = "black")+
guides(color = FALSE)+
ggtitle(bin2plot)+
theme(axis.text=element_text(size=14), axis.title=element_text(size=20),
title=element_text(size=20), legend.text=element_text(size=14),
legend.background = element_rect(fill="transparent"),
axis.text.x = element_text(angle = 45, hjust = 1),
strip.text.y=element_text(size=14))+
ylab("Density")+
xlab("Nucleotide identity (%)")+
xlim(75,100)
print(p_blast_sdisc)
}
# Plot one combined figure with proportions normalized for genome size
p_blast_sdisc_merged <- blast_df_sum %>%
ggplot(aes(x = bin_xcoord, ..density.., color = bin))+
theme_gray()+
scale_color_brewer("", palette = "Paired")+
facet_wrap(~Sample, nrow = 4)+
geom_density()+
ggtitle("Merged bins per sample")+
theme(axis.text=element_text(size=14), axis.title=element_text(size=20),
title=element_text(size=20), legend.text=element_text(size=14),
legend.background = element_rect(fill="transparent"),
axis.text.x = element_text(angle = 45, hjust = 1),
strip.text.y=element_text(size=14),
legend.position = "bottom")+
ylab("Reads per Mbp")+
xlab("Nucleotide identity (%)")+
xlim(75,100)+
guides(color=guide_legend(ncol=3))
print(p_blast_sdisc_merged)
# Plot for most abundant bin (B63)
p_blast_sdisc_B63 <- blast_df_sum %>%
dplyr::filter(bin == "B63_Su13.BD.MM110.DCMD_rebin1") %>%
ggplot(aes(x = bin_xcoord, ..density..))+
theme_bw()+
# facet_grid(season~Site)+
geom_density(alpha = 0.4, size = 0.75, color = "#333333")+
scale_shape_manual("", values = c(21,22,24))+
scale_color_brewer("", palette = "Accent")+
guides(fill = FALSE)+
theme(axis.text.y=element_text(size=14), axis.title=element_text(size=20),
title=element_text(size=20), legend.text=element_text(size=14),
legend.background = element_rect(fill="transparent"),
axis.text.x = element_text(size = 14),
strip.text=element_text(size=16),
panel.grid.minor = element_blank(),
legend.position = "bottom")+
ylab("Density")+
xlab("Nucleotide identity (%)")+
xlim(75,100)
p_blast_sdisc_B63
blast_df_sum$new_bin_name <- factor(blast_df_sum$new_bin_name, levels =
c("MAG5.SP-M110-DD","MAG2.FA-MLB-SN",
"MAG3.FA-MLB-SN", "MAG4.FA-M110-DN",
"MAG1.FA-MLB-DN", "MAG10.SU-M15-SN",
"MAG6.SP-M15-SD", "MAG8.SU-M110-DCMD",
"MAG7.SU-MLB-SD","MAG9.SU-M15-SN"))
# Plot for most abundant bin (B63)
p_blast_all_violin <- blast_df_sum %>%
ggplot(aes(x = new_bin_name, y= Identity, fill = new_bin_name))+
theme_bw()+
# facet_grid(season~Site)+
# geom_density(alpha = 0.4, size = 0.75, color = "#333333")+
geom_violin(alpha = 0.4, size = 0.75, color = "#333333",
scale = "area", fill = "#333333")+
scale_color_brewer("", palette = "Paired")+
guides(fill = FALSE)+
# facet_grid(new_bin_name~., scales = "free")+
theme(axis.text.y=element_text(size=14), axis.title=element_text(size=20),
title=element_text(size=20), legend.text=element_text(size=14),
legend.background = element_rect(fill="transparent"),
axis.text.x = element_text(size = 14),
strip.text=element_text(size=16),
panel.grid.minor = element_blank(),
legend.position = "bottom")+
# scale_x_log10()+
ylab("% Identity")+
xlab("")+
scale_y_continuous(limits = c(90,100))+
guides(shape = FALSE)+
# stat_summary(fun.data=mean_sdl, fun.args = list(mult = 1),
# geom="pointrange", color="#333333", size = 1)+
coord_flip()
# xlim(90,100)
p_blast_all_violin
# Plot for most abundant bin (B63)
p_blast_all_dens <- blast_df_sum %>%
ggplot(aes(x = bin_xcoord, shape = Sample))+
theme_bw()+
# facet_grid(season~Site)+
geom_density(alpha = 0.4, size = 0.75, color = "#333333")+
# geom_violin(alpha = 0.4, size = 0.75, color = "#333333",
# scale = "area", fill = "#333333")+
scale_color_brewer("", palette = "Paired")+
guides(fill = FALSE)+
facet_grid(new_bin_name~., scales = "free")+
theme(axis.text.y=element_text(size=14), axis.title=element_text(size=20),
title=element_text(size=20), legend.text=element_text(size=14),
legend.background = element_rect(fill="transparent"),
axis.text.x = element_text(size = 14),
strip.text=element_text(size=16),
panel.grid.minor = element_blank(),
legend.position = "bottom")+
# scale_x_log10()+
ylab("% Identity")+
xlab("")+
scale_x_continuous(limits = c(90,100))+
guides(shape = FALSE)
# stat_summary(fun.data=mean_sdl, fun.args = list(mult = 1),
# geom="pointrange", color="#333333", size = 1)+
# coord_flip()
# xlim(90,100)
p_blast_all_dens
# Plot for most Muskegon Lake for all bins
p_blast_sdisc_ML <- blast_df_sum %>% dplyr::filter(Site == "Muskegon Lake") %>%
ggplot(aes(x = bin_xcoord, ..scaled.., fill = new_bin_name, group = Sample,
shape = Depth))+
theme_bw()+
facet_grid(new_bin_name~season, scales = "free_y")+
geom_density(alpha = 0.3, color = "black", size = 1)+
# geom_point(size = 3, alpha = 0.6)+
scale_shape_manual("", values = c(21,22,24))+
scale_fill_manual(values = fill_palette)+
guides(color = FALSE, fill = FALSE)+
# ggtitle(bin2plot)+
theme(axis.text=element_text(size=14), axis.title=element_text(size=20),
title=element_text(size=20), legend.text=element_text(size=14),
legend.background = element_rect(fill="transparent"),
axis.text.x = element_text(angle = 45, hjust = 1),
strip.text=element_text(size=14),
panel.grid.minor = element_blank(),
legend.position = "bottom")+
ylab("Scaled density")+
xlab("Nucleotide identity (%)")+
xlim(75,100)
p_blast_sdisc_ML
blast_df_sum_comp <- blast_df_sum %>% group_by(Sample, bin) %>% dplyr::count(bin_xcoord)
blast_df_sum_comp <- dplyr::left_join(blast_df_sum_comp, bin_size, by = c("bin" = "bins"))
blast_df_sum_comp <- left_join(blast_df_sum_comp, new_bin_names, by = c("bin" = "bins"))
blast_df_sum_comp$new_bin_name <- as.character(blast_df_sum_comp$new_bin_name)
blast_df_sum_comp$new_bin_name <- factor(blast_df_sum_comp$new_bin_name, levels =
c("MAG1.FA-MLB-DN","MAG2.FA-MLB-SN",
"MAG3.FA-MLB-SN", "MAG4.FA-M110-DN",
"MAG5.SP-M110-DD","MAG6.SP-M15-SD",
"MAG7.SU-MLB-SD","MAG8.SU-M110-DCMD",
"MAG9.SU-M15-SN","MAG10.SU-M15-SN"))
blast_df_sum_comp <- dplyr::left_join(blast_df_sum_comp, meta_blast, by = c("Sample" = "Sample_ID"))
# Reorder site factor
blast_df_sum_comp$Site <- as.character(blast_df_sum_comp$Site)
blast_df_sum_comp$Site <- gsub("Buoy","Muskegon Lake", blast_df_sum_comp$Site)
blast_df_sum_comp$Site <- gsub("110","Lake Michigan\nsite M110", blast_df_sum_comp$Site)
blast_df_sum_comp$Site <- gsub("15","Lake Michigan\nsite M15", blast_df_sum_comp$Site)
blast_df_sum_comp$Site <- factor(blast_df_sum_comp$Site,
levels = c("Muskegon Lake",
"Lake Michigan\nsite M15",
"Lake Michigan\nsite M110"))
blast_df_sum_comp$Depth <- as.character(blast_df_sum_comp$Depth)
blast_df_sum_comp$Depth <- factor(blast_df_sum_comp$Depth,
levels = c("Surface", "Mid", "Deep"))
blast_df_sum_comp$Season <- as.character(blast_df_sum_comp$Season)
blast_df_sum_comp$Season <- factor(blast_df_sum_comp$Season, levels = c("Spring", "Summer", "Fall"))
# Normalize mapped reads per sample based on genome size
blast_df_sum_comp <- blast_df_sum_comp %>% group_by(bin) %>%
mutate(n_norm = 1e6*n/bin_size)
blast_df_sum_comp <- left_join(blast_df_sum_comp, total_reads, by = c("Sample" = "sample"))
# Divide normalized reads by 1M (fixed blast census)
blast_df_sum_comp <- blast_df_sum_comp %>% mutate(n_norm_perc = 100*n_norm/1e6)# Plot % reads corrected for genome size over threshold of 0.95
id_thresh <- 95-0.25
map_disc_cum <- blast_df_sum_comp %>% dplyr::filter(bin_xcoord > id_thresh) %>% group_by(Sample, bin) %>%
mutate(cum_rel_reads_mapped = cumsum(n_norm_perc))%>%
dplyr::filter(bin_xcoord == 100-0.25)
sum_cum <- map_disc_cum %>% group_by(Sample, bin) %>% mutate(cum_bins_rel_reads_mapped = sum(cum_rel_reads_mapped))
p_sdisc_cum3 <- ggplot(map_disc_cum, aes(x = new_bin_name,
y = cum_rel_reads_mapped,
fill = new_bin_name))+
theme_bw()+
scale_fill_manual("", values = fill_palette)+
geom_jitter(size = 4, shape = 21, color = "black", alpha = 0.7, width = 0.15)+
theme(axis.text=element_text(size=14), axis.title=element_text(size=20),
title=element_text(size=20), legend.text=element_text(size=12),
legend.background = element_rect(fill="transparent"),
axis.text.x = element_blank(),
strip.text=element_text(size=14), legend.position = "bottom",
strip.background = element_rect(fill = adjustcolor("gray", 0.15)))+
ylab(paste0("Norm. relative abundance ( > ", id_thresh, "% NI)"))+
xlab("")+
guides(fill=guide_legend(nrow = 3))+
facet_grid(Season~Site, scales ="free")+
scale_y_continuous(labels=scaleFUN, limits = c(0,3))+
coord_trans(y = "sqrt")
print(p_sdisc_cum3)
# Plot % reads over threshold of 0.99
id_thresh <- 99-0.25
map_disc_cum <- blast_df_sum_comp %>% dplyr::filter(bin_xcoord > id_thresh) %>% group_by(Sample, bin) %>%
mutate(cum_rel_reads_mapped = cumsum(n_norm_perc))%>%
dplyr::filter(bin_xcoord == 100-0.25)
sum_cum <- map_disc_cum %>% group_by(Sample, bin) %>% mutate(cum_bins_rel_reads_mapped = sum(cum_rel_reads_mapped))
p_sdisc_cum4 <- ggplot(map_disc_cum, aes(x = new_bin_name,
y = cum_rel_reads_mapped,
fill = new_bin_name))+
theme_bw()+
scale_fill_manual("", values = fill_palette)+
geom_jitter(size = 4, shape = 21, color = "black", alpha = 0.7, width = 0.15)+
theme(axis.text=element_text(size=14), axis.title=element_text(size=20),
title=element_text(size=20), legend.text=element_text(size=12),
legend.background = element_rect(fill="transparent"),
axis.text.x = element_blank(),
strip.text=element_text(size=14), legend.position = "bottom",
strip.background = element_rect(fill = adjustcolor("gray", 0.15)))+
ylab(paste0("Norm. relative abundance ( > ", id_thresh, "% NI)"))+
xlab("")+
guides(fill=guide_legend(nrow = 3))+
facet_grid(Season~Site, scales ="free")+
scale_y_continuous(labels=scaleFUN, limits = c(0,1.25))+
coord_trans(y = "sqrt")
print(p_sdisc_cum4)
Compare blast to bwa results
# Compare blast inferred abundances with bwa inferred abundances
grid_arrange_shared_legend(p_abs2, p_sdisc_cum4, ncol = 2)
# Merge those two dataframes and make scatter plot highlighting the differences
df_map_merged_sm <- df_map_merged[, c("Bin","Sample","rel_norm_abundance", "Season", "Site")]
df_map_merged_sm$Sample <- gsub(".", "_", df_map_merged_sm$Sample, fixed = TRUE)
df_map_merged_sm <- data.frame(df_map_merged_sm,
interaction = interaction(df_map_merged_sm$Bin,
df_map_merged_sm$Sample),
method = "bwa")
map_disc_cum_sm <- map_disc_cum[, c("bin","Sample","cum_rel_reads_mapped")]
map_disc_cum_sm <- data.frame(map_disc_cum_sm,
interaction = interaction(map_disc_cum_sm$bin,
map_disc_cum_sm$Sample),
method = "blastn")
abs_merged <- left_join(df_map_merged_sm, map_disc_cum_sm, by = "interaction")
abs_merged <- abs_merged[,c("Bin","Sample.x","rel_norm_abundance","cum_rel_reads_mapped",
"Season", "Site")]
colnames(abs_merged) <- c("Bin","Sample",
"bwa_norm_abundance","blastn_norm_abundance",
"Season", "Site")
# Plot
p_scat_abund <- ggplot(abs_merged, aes(x = bwa_norm_abundance, y = blastn_norm_abundance))+
theme_bw()+
geom_point(size = 4, shape = 21, color = "black", alpha = 0.7,
aes(fill = Bin))+
scale_fill_manual("", values = fill_palette)+
theme(axis.text=element_text(size=14), axis.title=element_text(size=20),
title=element_text(size=20), legend.text=element_text(size=12),
legend.background = element_rect(fill="transparent"),
strip.text=element_text(size=14), legend.position = "bottom",
strip.background = element_rect(fill = adjustcolor("gray", 0.15)))+
ylab(paste0("Norm. relative abundance (blastn, > ", id_thresh, "% NI)"))+
xlab("Norm. relative abundance (BWA)")+
guides(fill=guide_legend(nrow = 3))+
facet_grid(Season~Site, scales ="free")+
scale_y_continuous(labels=scaleFUN)+
scale_x_continuous(labels=scaleFUN)+
geom_smooth(method = "lm", color = "black")
# coord_trans(x = "atanh", y = "atanh")
print(p_scat_abund)
8. DESMAN
# Import data
results_desm <- read.csv("./DESMAN/MAG8_scg_10_9/Gamma_meanR.csv", header = TRUE)
colnames(results_desm)[1] <- "Samples"
colnames(results_desm) <- gsub("X", "Strain", colnames(results_desm))
colnames(results_desm) <- c("Samples", "Strain 1", "Strain 2", "Strain 3", "Strain 4",
"Strain 5")
results_desm$Samples <- gsub(".sorted.MAG8","",results_desm$Samples, fixed = TRUE)
# Import library prep data
library_prep <- read.table("./mapping_files/Libraries.tsv",header = TRUE)
# Wide to long format
results_desm_long <- tidyr::gather(results_desm, Strain, Freq,
2:(ncol(results_desm)), factor_key=TRUE)
# Merge with metadata
results_desm_long <- dplyr::left_join(results_desm_long, library_prep, by = "Samples")
results_desm_long$Samples <- gsub(".A", "", results_desm_long$Samples, fixed = TRUE)
results_desm_long$Samples <- gsub(".C", "", results_desm_long$Samples, fixed = TRUE)
results_desm_long <- dplyr::left_join(results_desm_long, meta_em, by = c("Samples" = "Sample_ID"))
results_desm_long$Site <- as.character(results_desm_long$Site)
results_desm_long$Site <- gsub("110", "Lake Michigan\nsite M110", results_desm_long$Site)
results_desm_long$Site <- gsub("15", "Lake Michigan\nsite M15", results_desm_long$Site)
results_desm_long$Site <- gsub("Buoy", "Muskegon Lake", results_desm_long$Site)
results_desm_long$Site <- factor(results_desm_long$Site, levels = c("Muskegon Lake",
"Lake Michigan\nsite M15",
"Lake Michigan\nsite M110"))
results_desm_long$Season <- as.character(results_desm_long$Season)
results_desm_long$Season <- factor(results_desm_long$Season, levels = c("Spring", "Summer","Fall"))
# Factor to numeric
results_desm_long$Temperature..C. <- as.numeric(as.character(results_desm_long$Temperature..C.))
results_desm_long$PAR <- as.numeric(as.character(results_desm_long$PAR))
results_desm_long$Chl.Lab..ug.L. <- as.numeric(as.character(results_desm_long$Chl.Lab..ug.L.))
results_desm_long$DO.Probe..mg.L. <- as.numeric(as.character(results_desm_long$DO.Probe..mg.L.))
results_desm_long$TP.ug.L <- as.numeric(as.character(results_desm_long$TP.ug.L))
results_desm_long$DOC.mg.L <- as.numeric(as.character(results_desm_long$DOC.mg.L))Check correlation with environmental variables
desm_p1b <- ggplot(results_desm_long, aes(x = Season, y = Freq,
fill = Strain,
col = Strain))+
geom_point(color = "black", alpha = 0.7, size = 3, shape = 21)+
# scale_shape_manual("", values = c(24,22,21))+
# scale_shape_manual(values = c(21,22,23))+
geom_boxplot(alpha = 0.3, col = "black", outlier.shape = NA, size = 0.3)+
scale_fill_brewer("", palette = "Accent")+
scale_color_brewer("",palette = "Accent")+
theme_bw()+
# geom_point(size = 4, color = "black", alpha = 0.7)+
# scale_shape_manual(values = c(21,24,23))+
# geom_boxplot(alpha = 0.4)+
theme(axis.text=element_text(size=12), axis.title=element_text(size=16),
title=element_text(size=16), legend.text=element_text(size=10),
legend.background = element_rect(fill="transparent"),
axis.text.x = element_text(size = 12, angle = 45, hjust = 1),
strip.text=element_text(size=12), legend.position = "bottom",
strip.background = element_rect(fill = adjustcolor("gray", 0.15)))+
ylab(paste0("Strain frequency"))+
guides(fill=FALSE)+
facet_grid(Site~Strain, scales ="free")+
xlab("")+
scale_y_continuous(labels=scaleFUN, limits = c(0,1))+
# geom_smooth(se = FALSE)+
guides(fill=FALSE)
# coord_trans(y = "sqrt")
print(desm_p1b)
# Just depth profile for M110 station
desm_p2b <- results_desm_long %>% dplyr::filter(Site == "Lake Michigan\nsite M110"&
Depth != "Mid") %>%
ggplot(., aes(x = Season, y = Freq,
fill = Strain,
col = Strain))+
geom_point(color = "black", alpha = 0.7, size = 4, shape = 21)+
scale_fill_brewer("", palette = "Accent")+
scale_color_brewer("",palette = "Accent")+
theme_bw()+
# geom_point(size = 4, color = "black", alpha = 0.7)+
# scale_shape_manual(values = c(21,24,23))+
# geom_boxplot(alpha = 0.4)+
theme(axis.text=element_text(size=12), axis.title=element_text(size=16),
title=element_text(size=16), legend.text=element_text(size=10),
legend.background = element_rect(fill="transparent"),
axis.text.x = element_text(size = 12, angle = 45, hjust = 1),
strip.text=element_text(size=12), legend.position = "bottom",
strip.background = element_rect(fill = adjustcolor("gray", 0.15)))+
ylab(paste0("Strain frequency"))+
guides(fill=FALSE)+
facet_grid(Depth~Strain, scales ="free")+
xlab("")+
scale_y_continuous(labels=scaleFUN, limits = c(0,1))+
# geom_smooth(se = FALSE)+
guides(fill=FALSE)
# coord_trans(y = "sqrt")
print(desm_p2b)
# Make correlation plot of strain variances alone
p.mat_t <- cor.mtest(results_desm[-1],method = "pearson", use="pairwise")
corrplot(cor(results_desm[-1], method = "pearson"), order = "hclust", addrect = 2,
tl.col="black", tl.srt=45,
sig.level = 0.05,
diag=TRUE)
print(cor(results_desm[-1], method = "pearson"))## Strain 1 Strain 2 Strain 3 Strain 4 Strain 5
## Strain 1 1.0000000 -0.8202727 -0.41160714 0.8601689 -0.61595753
## Strain 2 -0.8202727 1.0000000 0.10207177 -0.8626241 0.48346894
## Strain 3 -0.4116071 0.1020718 1.00000000 -0.2628128 -0.06622574
## Strain 4 0.8601689 -0.8626241 -0.26281282 1.0000000 -0.76358889
## Strain 5 -0.6159575 0.4834689 -0.06622574 -0.7635889 1.00000000results_desm$Samples <- gsub(".A", "", results_desm$Samples, fixed = TRUE)
results_desm$Samples <- gsub(".C", "", results_desm$Samples, fixed = TRUE)
# Make correlation plot of strain variances with environmental parameters
results_desm_env <- left_join(results_desm, meta_em, by = c("Samples" = "Sample_ID"))
results_desm_env <- results_desm_env %>%
dplyr::select(Samples:`Strain 5`, Temperature..C., `Cond..µS.cm.`,
TP.ug.L, DOC.mg.L, PAR, DO.Probe..mg.L.,
Chl.Lab..ug.L.) %>%
apply(., 2, function(x) as.numeric(as.character(x))) %>% data.frame()
corrplot(cor(results_desm_env[-1], method = "pearson", use="pairwise"),
addrect = 2, order = "hclust",
tl.col="black", tl.srt=45,
sig.level = 0.05,
diag=TRUE)
print(cor(results_desm_env[-1], method = "pearson", use="pairwise")[1:5, ])## Strain.1 Strain.2 Strain.3 Strain.4 Strain.5
## Strain.1 1.0000000 -0.8202727 -0.41160715 0.8601689 -0.61595750
## Strain.2 -0.8202727 1.0000000 0.10207178 -0.8626241 0.48346889
## Strain.3 -0.4116071 0.1020718 1.00000000 -0.2628128 -0.06622572
## Strain.4 0.8601689 -0.8626241 -0.26281282 1.0000000 -0.76358887
## Strain.5 -0.6159575 0.4834689 -0.06622572 -0.7635889 1.00000000
## Temperature..C. Cond..µS.cm. TP.ug.L DOC.mg.L PAR
## Strain.1 0.34481960 -0.4860871 -0.3861337 0.13104041 0.14442728
## Strain.2 -0.07366979 0.5616745 0.6354257 0.13681522 -0.18508745
## Strain.3 -0.90955092 -0.3110741 0.1936932 -0.07672245 -0.36314128
## Strain.4 0.09272609 -0.3154406 -0.5826171 -0.31330645 -0.04343812
## Strain.5 0.25378964 0.4672220 0.7396795 0.30023012 0.20455637
## DO.Probe..mg.L. Chl.Lab..ug.L.
## Strain.1 -0.03438865 -0.25151040
## Strain.2 -0.06078842 0.38969012
## Strain.3 0.42427714 0.04255277
## Strain.4 0.21315075 -0.40976670
## Strain.5 -0.50569631 0.65365968colnames(results_desm_env)[-1]## [1] "Strain.1" "Strain.2" "Strain.3"
## [4] "Strain.4" "Strain.5" "Temperature..C."
## [7] "Cond..µS.cm." "TP.ug.L" "DOC.mg.L"
## [10] "PAR" "DO.Probe..mg.L." "Chl.Lab..ug.L."print(cor.mtest(results_desm_env[-1],method = "pearson", use="pairwise")$p)## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.000000e+00 9.239903e-07 4.567366e-02 7.108310e-08 1.352611e-03
## [2,] 9.239903e-07 0.000000e+00 6.350749e-01 5.920205e-08 1.669098e-02
## [3,] 4.567366e-02 6.350749e-01 0.000000e+00 2.147037e-01 7.584960e-01
## [4,] 7.108310e-08 5.920205e-08 2.147037e-01 0.000000e+00 1.417411e-05
## [5,] 1.352611e-03 1.669098e-02 7.584960e-01 1.417411e-05 0.000000e+00
## [6,] 1.160468e-01 7.445643e-01 4.534080e-09 6.814959e-01 2.544210e-01
## [7,] 7.801333e-02 3.659950e-02 2.790026e-01 2.719500e-01 9.209330e-02
## [8,] 1.551432e-01 1.090805e-02 4.891349e-01 2.265657e-02 1.621215e-03
## [9,] 6.415704e-01 6.268178e-01 7.858018e-01 2.555027e-01 2.769303e-01
## [10,] 5.935779e-01 4.925486e-01 1.668248e-01 8.730920e-01 4.472885e-01
## [11,] 8.823569e-01 7.935129e-01 5.524854e-02 3.535678e-01 1.935020e-02
## [12,] 3.473915e-01 1.356871e-01 8.756597e-01 1.149672e-01 6.027313e-03
## [,6] [,7] [,8] [,9] [,10]
## [1,] 1.160468e-01 0.07801333 1.551432e-01 0.64157041 0.59357792
## [2,] 7.445643e-01 0.03659950 1.090805e-02 0.62681784 0.49254864
## [3,] 4.534080e-09 0.27900255 4.891349e-01 0.78580179 0.16682479
## [4,] 6.814959e-01 0.27195000 2.265657e-02 0.25550268 0.87309199
## [5,] 2.544210e-01 0.09209330 1.621215e-03 0.27693035 0.44728845
## [6,] 0.000000e+00 0.46841256 3.470042e-01 0.42301090 0.02338948
## [7,] 4.684126e-01 0.00000000 6.918041e-01 0.97959502 0.69252132
## [8,] 3.470042e-01 0.69180406 0.000000e+00 0.29030153 0.72210624
## [9,] 4.230109e-01 0.97959502 2.903015e-01 0.00000000 0.14215923
## [10,] 2.338948e-02 0.69252132 7.221062e-01 0.14215923 0.00000000
## [11,] 7.414131e-03 0.01218690 2.385871e-01 0.06807210 0.45911626
## [12,] 6.392885e-01 0.63632415 4.665197e-07 0.08948707 0.87398504
## [,11] [,12]
## [1,] 0.882356890 3.473915e-01
## [2,] 0.793512882 1.356871e-01
## [3,] 0.055248542 8.756597e-01
## [4,] 0.353567755 1.149672e-01
## [5,] 0.019350202 6.027313e-03
## [6,] 0.007414131 6.392885e-01
## [7,] 0.012186900 6.363242e-01
## [8,] 0.238587083 4.665197e-07
## [9,] 0.068072096 8.948707e-02
## [10,] 0.459116263 8.739850e-01
## [11,] 0.000000000 6.159640e-01
## [12,] 0.615964034 0.000000e+00desm_p2 <- ggplot(results_desm_long, aes(x = Temperature..C., y = Freq,
fill = Strain,
col = Strain))+
geom_point(color = "black", alpha = 0.7, shape = 21, size = 4, col = "black")+
scale_fill_brewer("", palette = "Accent")+
scale_color_brewer("", palette = "Accent")+
theme_bw()+
# geom_point(size = 4, color = "black", alpha = 0.7)+
# scale_shape_manual(values = c(21,24,23))+
# geom_boxplot(alpha = 0.4)+
theme(axis.text=element_text(size=12), axis.title=element_text(size=16),
title=element_text(size=16), legend.text=element_text(size=10),
legend.background = element_rect(fill="transparent"),
axis.text.x = element_text(size = 12),
strip.text=element_text(size=12), legend.position = "bottom",
strip.background = element_rect(fill = adjustcolor("gray", 0.15)))+
ylab(paste0("Strain frequency"))+
guides(fill=FALSE)+
# facet_grid(Season~Site, scales ="free")+
xlab("Temperature (°C)")+
scale_y_continuous(labels=scaleFUN, limits = c(0,1))+
geom_smooth(se = FALSE)+
guides(fill=guide_legend(nrow=2))
# coord_trans(y = "sqrt")
desm_p3 <- results_desm_long %>% mutate(Chl.Lab..ug.L. =
as.numeric(as.character(Chl.Lab..ug.L.))) %>%
ggplot(., aes(x = Chl.Lab..ug.L. , y = Freq,
fill = Strain,
col = Strain))+
geom_point(color = "black", alpha = 0.7, shape = 21, size = 4, col = "black")+
scale_fill_brewer("", palette = "Accent")+
scale_color_brewer("", palette = "Accent")+
theme_bw()+
# geom_point(size = 4, color = "black", alpha = 0.7)+
# scale_shape_manual(values = c(21,24,23))+
# geom_boxplot(alpha = 0.4)+
theme(axis.text=element_text(size=12), axis.title=element_text(size=16),
title=element_text(size=16), legend.text=element_text(size=10),
legend.background = element_rect(fill="transparent"),
axis.text.x = element_text(size = 12),
strip.text=element_text(size=12), legend.position = "bottom",
strip.background = element_rect(fill = adjustcolor("gray", 0.15)))+
ylab(paste0("Strain frequency"))+
guides(fill=FALSE)+
# facet_grid(Season~Site, scales ="free")+
xlab("Chl a concentration (µg/L)")+
scale_y_continuous(labels=scaleFUN, limits = c(0,1))+
geom_smooth(se = FALSE, span = 1)+
guides(fill=guide_legend(nrow=2))
# coord_trans(y = "sqrt")
desm_p4 <- results_desm_long %>% mutate(TDP.SRP = as.numeric(as.character(TDP.SRP))) %>%
ggplot(., aes(x = TDP.SRP, y = Freq,
fill = Strain,
col = Strain))+
geom_point(color = "black", alpha = 0.7, shape = 21, size = 4, col = "black")+
scale_fill_brewer("", palette = "Accent")+
scale_color_brewer("",palette = "Accent")+
theme_bw()+
# geom_point(size = 4, color = "black", alpha = 0.7)+
# scale_shape_manual(values = c(21,24,23))+
# geom_boxplot(alpha = 0.4)+
theme(axis.text=element_text(size=12), axis.title=element_text(size=16),
title=element_text(size=16), legend.text=element_text(size=10),
legend.background = element_rect(fill="transparent"),
axis.text.x = element_text(size = 12),
strip.text=element_text(size=12), legend.position = "bottom",
strip.background = element_rect(fill = adjustcolor("gray", 0.15)))+
ylab(paste0("Strain frequency"))+
guides(fill=FALSE)+
# facet_grid(Season~Site, scales ="free")+
xlab("Soluble Reactive Phosphorus (SRP - µg/L)")+
scale_y_continuous(labels=scaleFUN, limits = c(0,1))+
geom_smooth(se = FALSE, span = 0.75)+
guides(fill=guide_legend(nrow=2))
# coord_trans(y = "sqrt")
desm_p5 <- results_desm_long %>% mutate(PON.mg.L = as.numeric(as.character(PON.mg.L))) %>%
ggplot(., aes(x = PON.mg.L, y = Freq,
fill = Strain,
col = Strain))+
geom_point(color = "black", alpha = 0.7, shape = 21, size = 4, col = "black")+
scale_fill_brewer("", palette = "Accent")+
scale_color_brewer("",palette = "Accent")+
theme_bw()+
# geom_point(size = 4, color = "black", alpha = 0.7)+
# scale_shape_manual(values = c(21,24,23))+
# geom_boxplot(alpha = 0.4)+
theme(axis.text=element_text(size=12), axis.title=element_text(size=16),
title=element_text(size=16), legend.text=element_text(size=10),
legend.background = element_rect(fill="transparent"),
axis.text.x = element_text(size = 12),
strip.text=element_text(size=12), legend.position = "bottom",
strip.background = element_rect(fill = adjustcolor("gray", 0.15)))+
ylab(paste0("Strain frequency"))+
guides(fill=FALSE)+
# facet_grid(Season~Site, scales ="free")+
xlab("Particulate Organic Nitrogen (PON - mg/L)")+
scale_y_continuous(labels=scaleFUN, limits = c(0,1))+
geom_smooth(se = FALSE, span = 0.75)+
guides(fill=guide_legend(nrow=2))
# coord_trans(y = "sqrt")
desm_p6 <- results_desm_long %>% mutate(POC.mg.L = as.numeric(as.character(POC.mg.L))) %>%
ggplot(., aes(x = POC.mg.L, y = Freq,
fill = Strain,
col = Strain))+
geom_point(color = "black", alpha = 0.7, shape = 21, size = 4, col = "black")+
scale_fill_brewer("", palette = "Accent")+
scale_color_brewer("",palette = "Accent")+
theme_bw()+
# geom_point(size = 4, color = "black", alpha = 0.7)+
# scale_shape_manual(values = c(21,24,23))+
# geom_boxplot(alpha = 0.4)+
theme(axis.text=element_text(size=12), axis.title=element_text(size=16),
title=element_text(size=16), legend.text=element_text(size=10),
legend.background = element_rect(fill="transparent"),
axis.text.x = element_text(size = 12),
strip.text=element_text(size=12), legend.position = "bottom",
strip.background = element_rect(fill = adjustcolor("gray", 0.15)))+
ylab(paste0("Strain frequency"))+
guides(fill=FALSE)+
# facet_grid(Season~Site, scales ="free")+
xlab("Particulate Organic Carbon (POC - mg/L)")+
scale_y_continuous(labels=scaleFUN, limits = c(0,1))+
geom_smooth(se = FALSE, span = 0.75)+
guides(fill=guide_legend(nrow=2))
# coord_trans(y = "sqrt")
# Combine all exploratory plots
cowplot::plot_grid(desm_p3, desm_p4, desm_p5, desm_p6, ncol = 2, align = 'hv')## `geom_smooth()` using method = 'loess'
## `geom_smooth()` using method = 'loess'
## `geom_smooth()` using method = 'loess'
## `geom_smooth()` using method = 'loess'
Diversity plots
alpha_div_desman <- data.frame(Sample = results_desm[,1],
D2 = vegan::diversity(results_desm[,2:ncol(results_desm)], index = "invsimpson")
)
alpha_div_desman <- left_join(alpha_div_desman, meta_em,
by = c("Sample" = "Sample_ID"))
alpha_div_desman$Site <- as.character(alpha_div_desman$Site)
alpha_div_desman$Site <- gsub("110", "Lake Michigan\nsite M110", alpha_div_desman$Site)
alpha_div_desman$Site <- gsub("15", "Lake Michigan\nsite M15", alpha_div_desman$Site)
alpha_div_desman$Site <- gsub("Buoy", "Muskegon Lake", alpha_div_desman$Site)
alpha_div_desman$Site <- factor(alpha_div_desman$Site, levels = c("Muskegon Lake",
"Lake Michigan\nsite M15",
"Lake Michigan\nsite M110"))
alpha_div_desman$Season <- as.character(alpha_div_desman$Season)
alpha_div_desman$Season <- factor(alpha_div_desman$Season, levels = c("Spring", "Summer","Fall"))
p_alpha_desman <- alpha_div_desman %>%
ggplot(., aes(x = Season, y = D2, fill = Season, shape = Depth))+
geom_point(size = 4) +
scale_shape_manual("", values = c(21,22,24))+
facet_grid(.~Site)+
theme_bw()+
scale_fill_brewer(palette = "Accent")+
labs(y = expression("LimB strain diversity (D"[2]*")"), x = "")+
theme(axis.text=element_text(size=12), axis.title=element_text(size=16),
title=element_text(size=16), legend.text=element_text(size=16),
strip.text = element_text(size = 16))+
ylim(0,4.2)+
guides(fill = guide_legend(title = "", override.aes = list(size = 4, shape = 21),
nrow = 4)
)
print(p_alpha_desman)
# Make PCoA
variants <- results_desm[,2:ncol(results_desm)]
rownames(variants) <- results_desm[, 1]
taxa <- as.matrix(data.frame(variant_name = colnames(results_desm)[2:ncol(results_desm)]))
rownames(taxa) <- colnames(results_desm)[2:ncol(results_desm)]
# Store as phyloseq object
physeq_desm <- phyloseq(otu_table(variants, taxa_are_rows = FALSE),
tax_table(taxa)
)
# Run beta diversity analysis on 16s data
pcoa <- ordinate(
physeq = physeq_desm,
method = "PCoA",
distance = "bray",
correction = "lingoes",
k=2
)
pcoa.df <- data.frame(Samples = sample_names(physeq_desm), pcoa$vectors)
var <- round(pcoa$values$Eigenvalues/sum(pcoa$values$Eigenvalues)*100,1)
# Add metadata
pcoa.df$Samples <- gsub(".A", "", pcoa.df$Samples, fixed = TRUE)
pcoa.df$Samples <- gsub(".C", "", pcoa.df$Samples, fixed = TRUE)
pcoa.df <- dplyr::left_join(pcoa.df, meta_em, by = c("Samples" = "Sample_ID"))
pcoa.df$Site <- as.character(pcoa.df$Site)
pcoa.df$Site <- gsub("110", "Lake Michigan - M110", pcoa.df$Site)
pcoa.df$Site <- gsub("15", "Lake Michigan - M15", pcoa.df$Site)
pcoa.df$Site <- gsub("Buoy", "Muskegon Lake", pcoa.df$Site)
pcoa.df$Site <- factor(pcoa.df$Site, levels = c("Muskegon Lake",
"Lake Michigan - M15",
"Lake Michigan - M110"))
pcoa.df$Season <- as.character(pcoa.df$Season)
pcoa.df$Season <- factor(pcoa.df$Season, levels = c("Spring", "Summer","Fall"))
# Run permanova
x <- results_desm[,2:ncol(results_desm)]
rownames(x) <- results_desm[, 1]
dist_desm <- vegan::vegdist(x)
perm <- permute::how(nperm = 999)
permanova_desm <- vegan::adonis(dist_desm ~ Season * Site,
data = meta_em,
permutations = perm)
# Plot PCoA
my_grob = grobTree(textGrob(bquote(paste(r[Season]^2 ==
.(round(100 * permanova_desm$aov.tab[1, 5], 1)),
"%")), x = 0.7, y = 0.95, hjust = 0, gp = gpar(col = "black",
fontsize = 14, fontface = "italic")))
my_grob2 = grobTree(textGrob(bquote(paste(r[Site]^2 ==
.(format(round(100 * permanova_desm$aov.tab[2, 5],
1), nsmall = 1)), "%")), x = 0.7, y = 0.87,
hjust = 0, gp = gpar(col = "black", fontsize = 14,
fontface = "italic")))
my_grob3 = grobTree(textGrob(bquote(paste(r[Season:Site]^2 ==
.(round(100 * permanova_desm$aov.tab[3, 5], 1)),
"%")), x = 0.7, y = 0.79, hjust = 0, gp = gpar(col = "black",
fontsize = 14, fontface = "italic")))
pcoa.ord <- ggplot(data=pcoa.df, aes(x=Axis.1, y=Axis.2, shape = Site))+
scale_shape_manual(values=c(21,23, 24))+
geom_point(size=7, alpha=0.7, aes(fill = Season))+
theme_bw()+
scale_fill_brewer(palette = "Accent")+
labs(x = paste0("PCoA axis 1 (",var[1], "%)"),
y = paste0("PCoA axis 2 (",var[2], "%)"), fill="", shape = "", colour="")+
theme(axis.text=element_text(size=14), axis.title=element_text(size=16),
title=element_text(size=16), legend.text=element_text(size=16))+
guides(fill = guide_legend(override.aes = list(shape = 22)))+
annotation_custom(my_grob)+
annotation_custom(my_grob2)+
annotation_custom(my_grob3)+
ggtitle("PCoA of resolved MAG8 variant composition")
print(pcoa.ord)
Gene assignment
# Import etaS.csv which contains gene assignment
geneAssign_df <- read.csv("./DESMAN/MAG8_scg_10_9/MAG8.SU-M110-DCMDetaS_df.csv",
header = TRUE)
colnames(geneAssign_df) <- gsub("X", "Variant", colnames(geneAssign_df))
geneAssign_df <- geneAssign_df[, c(1,6,8,9,10,11)]
colnames(geneAssign_df) <- c("Gene", "Strain 1", "Strain 2", "Strain 3", "Strain 4",
"Strain 5")
# Check for presence of core gene presence in strains
core_cogs_lim <- unique(read.csv("./DESMAN/MAG8_scg_10_9/Filtered_Tau_star.csv")[,1])
# Check predicted presence of core genes
geneAssign_df %>% dplyr::filter(Gene %in% core_cogs_lim)## Gene Strain 1 Strain 2 Strain 3 Strain 4 Strain 5
## 1 contig_20975_3 1 1 1 0 1
## 2 contig_20975_2 1 1 1 1 1
## 3 contig_20975_5 1 1 1 1 1
## 4 contig_20975_4 1 1 1 1 1
## 5 contig_20975_7 1 1 1 1 1
## 6 contig_21188_3 1 1 1 1 1
## 7 contig_21188_1 1 1 1 1 1
## 8 contig_21188_4 1 1 1 1 1
## 9 contig_21013_3 1 1 1 0 1
## 10 contig_21099_3 1 1 1 1 1
## 11 contig_20877_2 1 1 1 1 1
## 12 contig_20877_3 1 1 1 1 1
## 13 contig_20984_3 1 1 1 0 1
## 14 contig_20984_7 1 1 1 1 1
## 15 contig_20984_6 1 1 1 1 1
## 16 contig_20984_8 1 1 1 1 1
## 17 contig_20881_1 1 1 1 1 1
## 18 contig_20915_2 1 1 1 0 1# Wide to long format
geneAssign_df <- tidyr::gather(geneAssign_df, Variant, Presence,
`Strain 1`:`Strain 5`, factor_key=TRUE)
# Import blast results
blast_desman <- read.delim("./DESMAN/MAG8_scg_10_9/desman2img_genes.blast",
header = FALSE)
colnames(blast_desman) <- c("qseqid", "sseqid", "pident", "length", "mismatch",
"gapopen", "qstart", "qend", "sstart", "send",
"evalue", "bitscore")
blast_desman_sb <- blast_desman %>% dplyr::select("qseqid", "sseqid") %>% distinct()
# Merge IMG gene IDs with
geneAssign_df <- left_join(geneAssign_df, blast_desman_sb, by = c("Gene" = "qseqid"))
geneAssign_df$sseqid <- as.character(geneAssign_df$sseqid)
# Merge with functional annotation (KO & COG)
geneAssign_df_annot <- left_join(geneAssign_df, merged_file_annot,
by = c("sseqid" = "gene_oid"))
# Merge with metaT data
geneAssign_df_annot <- left_join(geneAssign_df, expr_cov_long,
by = c("sseqid" = "gene_oid"))Control core gene presence in strains
upset diagram
# Long to wide
geneAssign_df_wide <- tidyr::spread(geneAssign_df, Variant, Presence)
# Plot
upset(geneAssign_df_wide, sets = c("Strain 1", "Strain 2", "Strain 3", "Strain 4",
"Strain 5"), mb.ratio = c(0.6, 0.45),
order.by = "freq", point.size = 3.5,
mainbar.y.label = "Gene intersections", sets.x.label = "Number of genes",
text.scale = c(1.5, 1.5, 1.5, 1.4, 2, 1.5),
scale.intersections = "log2",
keep.order = FALSE,
line.size = 1.5)
upset(geneAssign_df_wide, sets = c("Strain 1", "Strain 4"), mb.ratio = c(0.6, 0.45),
order.by = "freq", point.size = 3.5,
mainbar.y.label = "Gene intersections", sets.x.label = "Number of genes",
text.scale = c(1.5, 1.5, 1.5, 1.4, 2, 2),
scale.intersections = "log2",
keep.order = FALSE,
line.size = 1.5)
upset(geneAssign_df_wide, sets = c("Strain 2", "Strain 5"), mb.ratio = c(0.6, 0.45),
order.by = "freq", point.size = 3.5,
mainbar.y.label = "Gene intersections", sets.x.label = "Number of genes",
text.scale = c(1.5, 1.5, 1.5, 1.4, 2, 2),
scale.intersections = "log2",
keep.order = FALSE,
line.size = 1.5)
# Indicate site and season-specific differential abundant genes
geneAssign_df_wide <- left_join(geneAssign_df_wide, res_deseq,
by = c("sseqid" = "gene_oid"))Site downregulation
# Make upset plot for site-specific downregulation controlled for seasons
geneAssign_df_wide %>% dplyr::filter(regulation == "downregulation",
Design == "~ Season + Site") %>%
upset(., sets = c("Strain 1", "Strain 2", "Strain 3", "Strain 4",
"Strain 5"), mb.ratio = c(0.55, 0.45),
order.by = "freq", number.angles = 30, point.size = 3.5,
mainbar.y.label = "Gene intersections", sets.x.label = "Number of genes",
text.scale = c(1.5, 1.5, 1.5, 1.4, 2, 0.75),
show.numbers = FALSE,
scale.intersections = "log2",
keep.order = FALSE,
boxplot.summary = "log2FoldChange",
line.size = NA)
Site upregulation
# Make upset plot for site-specific upregulation controlled for seasons
geneAssign_df_wide %>% dplyr::filter(regulation == "upregulation",
Design == "~ Season + Site") %>%
upset(., sets = c("Strain 1", "Strain 2", "Strain 3", "Strain 4",
"Strain 5"), mb.ratio = c(0.55, 0.45),
order.by = "freq", number.angles = 30, point.size = 3.5,
mainbar.y.label = "Gene intersections", sets.x.label = "Number of genes",
text.scale = c(1.5, 1.5, 1.5, 1.4, 2, 0.75),
show.numbers = FALSE,
scale.intersections = "log2",
keep.order = FALSE,
boxplot.summary = "log2FoldChange",
line.size = NA)
Seasonal downregulation
# Make upset plot for seasonal downregulation controlled for site
geneAssign_df_wide %>% dplyr::filter(regulation == "downregulation",
Design == "~ Site + Season") %>%
upset(., sets = c("Strain 1", "Strain 2", "Strain 3", "Strain 4",
"Strain 5"), mb.ratio = c(0.55, 0.45),
order.by = "freq", number.angles = 30, point.size = 3.5,
mainbar.y.label = "Gene intersections", sets.x.label = "Number of genes",
text.scale = c(1.5, 1.5, 1.5, 1.4, 2, 0.75),
show.numbers = FALSE,
scale.intersections = "log2",
keep.order = FALSE,
boxplot.summary = "log2FoldChange",
line.size = NA)
Depth downregulation
# Make upset plot for depth-specific downregulation controlled for seasons
geneAssign_df_wide %>% dplyr::filter(regulation == "downregulation",
Design == "~ Site + Depth") %>%
upset(., sets = c("Strain 1", "Strain 2", "Strain 3", "Strain 4",
"Strain 5"), mb.ratio = c(0.55, 0.45),
order.by = "freq", number.angles = 30, point.size = 3.5,
mainbar.y.label = "Gene intersections", sets.x.label = "Number of genes",
text.scale = c(1.5, 1.5, 1.5, 1.4, 2, 0.75),
show.numbers = FALSE,
scale.intersections = "log2",
keep.order = FALSE,
boxplot.summary = "log2FoldChange",
line.size = NA)
Depth upregulation
# Make upset plot for depth-specific upregulation controlled for sites
geneAssign_df_wide %>% dplyr::filter(regulation == "upregulation",
Design == "~ Site + Depth") %>%
upset(., sets = c("Strain 1", "Strain 2", "Strain 3", "Strain 4",
"Strain 5"), mb.ratio = c(0.55, 0.45),
order.by = "freq", number.angles = 30, point.size = 3.5,
mainbar.y.label = "Gene intersections", sets.x.label = "Number of genes",
text.scale = c(1.5, 1.5, 1.5, 1.4, 2, 0.75),
show.numbers = FALSE,
scale.intersections = "log2",
keep.order = FALSE,
boxplot.summary = "log2FoldChange",
line.size = NA)
Transcript assignment
# Extract upregulated genes of depth comparison
geneAssign_df_wide_transcripts_depth <- geneAssign_df_wide %>%
dplyr::filter(sseqid %in% unique(MAG8_deseq_results_depth_fin$gene_oid)) %>%
dplyr::select(sseqid:`Strain 5`)
geneAssign_df_wide_transcripts_depth <- MAG8_deseq_results_depth_fin %>%
dplyr::select(log2FoldChange, gene_oid, regulation) %>%
left_join(geneAssign_df_wide_transcripts_depth, ., by = c("sseqid" = "gene_oid"))
# Extract upregulated genes of surface comparison
geneAssign_df_wide_transcripts_temp <- geneAssign_df_wide %>%
dplyr::filter(sseqid %in% unique(MAG8_deseq_results_temp_fin$gene_oid)) %>%
dplyr::select(sseqid:`Strain 5`)
geneAssign_df_wide_transcripts_temp <- MAG8_deseq_results_temp_fin %>%
dplyr::select(log2FoldChange, gene_oid, regulation) %>%
left_join(geneAssign_df_wide_transcripts_temp, ., by = c("sseqid" = "gene_oid"))
# Evaluate distribution of transcripts expression over strains
# Plot upregulated genes
geneAssign_df_wide_transcripts_depth %>%
unique() %>%
upset(.,
sets = c("Strain 1", "Strain 2","Strain 3", "Strain 4", "Strain 5"),
mb.ratio = c(0.55, 0.45),
order.by = "freq", point.size = 3.5,
mainbar.y.label = "Gene intersections", sets.x.label = "Number of genes",
text.scale = c(1.5, 1.5, 1.5, 1.4, 2, 3),
scale.intersections = "log2",
keep.order = FALSE,
line.size = 1.2)
geneAssign_df_wide_transcripts_depth %>%
dplyr::filter(regulation == "Downregulated") %>%
unique() %>%
summary()## sseqid Strain 1 Strain 2 Strain 3
## Length:45 Min. :0.0000 Min. :0.0000 Min. :0.0000
## Class :character 1st Qu.:1.0000 1st Qu.:0.0000 1st Qu.:1.0000
## Mode :character Median :1.0000 Median :1.0000 Median :1.0000
## Mean :0.9333 Mean :0.5778 Mean :0.9778
## 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:1.0000
## Max. :1.0000 Max. :1.0000 Max. :1.0000
## Strain 4 Strain 5 log2FoldChange regulation
## Min. :0.0000 Min. :0.0000 Min. :-4.535 Upregulated : 0
## 1st Qu.:0.0000 1st Qu.:1.0000 1st Qu.:-3.143 Downregulated:45
## Median :1.0000 Median :1.0000 Median :-2.521
## Mean :0.6444 Mean :0.9111 Mean :-2.689
## 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:-2.179
## Max. :1.0000 Max. :1.0000 Max. :-1.684#
# geneAssign_df_wide_transcripts_temp %>% dplyr::filter(regulation == "Upregulated") %>%
# upset(.,
# sets = c("Strain 1", "Strain 2", "Strain 3", "Strain 4",
# "Strain 5"), mb.ratio = c(0.55, 0.45),
# order.by = "freq", number.angles = 30, point.size = 3.5,
# mainbar.y.label = "Gene intersections", sets.x.label = "Number of genes",
# text.scale = c(1.5, 1.5, 1.5, 1.4, 2, 0.75),
# show.numbers = FALSE,
# scale.intersections = "log2",
# keep.order = FALSE,
# line.size = 1.2,
# boxplot.summary = "log2FoldChange")
# geneAssign_df_wide_transcripts_temp %>% dplyr::filter(regulation == "Downregulated") %>%
# upset(.,
# sets = c("Strain 1", "Strain 2", "Strain 3", "Strain 4",
# "Strain 5"), mb.ratio = c(0.55, 0.45),
# order.by = "freq", number.angles = 30, point.size = 3.5,
# mainbar.y.label = "Gene intersections", sets.x.label = "Number of genes",
# text.scale = c(1.5, 1.5, 1.5, 1.4, 2, 0.75),
# show.numbers = FALSE,
# scale.intersections = "log2",
# keep.order = FALSE,
# line.size = 1.2,
# boxplot.summary = "log2FoldChange")Assign functions
# Assign interactins for each strain
geneAssign_df_wide_transcripts_depth$Strain_interaction <-
interaction(geneAssign_df_wide_transcripts_depth$`Strain 1`,
geneAssign_df_wide_transcripts_depth$`Strain 2`,
geneAssign_df_wide_transcripts_depth$`Strain 3`,
geneAssign_df_wide_transcripts_depth$`Strain 4`,
geneAssign_df_wide_transcripts_depth$`Strain 5`,
sep = "-")
# Assign functions
geneAssign_df_wide_transcripts_depth <- left_join(geneAssign_df_wide_transcripts_depth,
merged_file_annot,
by = c("sseqid" = "gene_oid"))
# Check for gene enrichment in interactions of strain 1 - 4 / strain 2 - 5 / strain 3
# these interactions correspond to interaction term 1-0-0-1-0 / 0-1-0-0-1 / 0-0-1-0-0
# geneAssign_df_wide_transcripts_depth %>%
# dplyr::filter(grepl("1-.-.-1-.", Strain_interaction))
geneAssign_df_wide_transcripts_depth %>%
# dplyr::filter(!grepl("0-0-0-0-0", Strain_interaction)) %>%
# dplyr::filter(grepl("1-0-1-0-0", Strain_interaction)) %>%
dplyr::select(sseqid:ko_name) %>%
# dplyr::filter(!is.na(ko_name)) %>%
unique() %>%
ggplot(., aes(x = Strain_interaction, y = log2FoldChange, fill = regulation))+
# geom_point(shape = 21, size = 4)+
geom_boxplot(alpha = 0.4, outlier.shape = NA)+
# geom_bar(stat = "count")+
scale_fill_brewer("", palette = "Set1")+
theme_bw()+
theme(axis.text.x=element_text(size=12, angle = 30, hjust = 1),
axis.text.y = element_text(size = 12),
axis.title=element_text(size=16),
title=element_text(size=16), legend.text=element_text(size=16),
strip.text = element_text(size = 16))+
guides(fill = FALSE)
# ylim(0,40)Black Queen Hypothesis?
geneAssign_df_annot_func <- geneAssign_df_wide[, 1:7]
geneAssign_df_annot_func$Strain_interaction <-
interaction(geneAssign_df_annot_func$`Strain 1`,
geneAssign_df_annot_func$`Strain 2`,
geneAssign_df_annot_func$`Strain 3`,
geneAssign_df_annot_func$`Strain 4`,
geneAssign_df_annot_func$`Strain 5`,
sep = "-")
# Annotate full accessory genomes
geneAssign_df_annot_func <- left_join(geneAssign_df_annot_func,
merged_file_annot, by = c("sseqid" = "gene_oid"))
# Compare only strain 1 and 4
p_unique_s1 <- geneAssign_df_annot_func %>%
dplyr::filter(grepl("1-.-.-0-.", Strain_interaction)) %>%
dplyr::select(sseqid, ko_name, ko_function_abbrev:ko_level_C) %>%
dplyr::filter(!is.na(ko_name)) %>%
unique() %>%
ggplot(., aes(x = ko_level_C, fill = ko_level_B))+
# geom_point()+
# geom_boxplot()+
geom_bar(stat = "count", col = "black")+
scale_fill_brewer("", palette = "Paired")+
theme_bw()+
theme(axis.text.x=element_text(size=12, angle = 30, hjust = 1),
axis.text.y = element_text(size = 12),
axis.title=element_text(size=16),
title=element_text(size=16), legend.text=element_text(size=16),
strip.text = element_text(size = 16))+
guides(fill = FALSE)
# ylim(0,40)
p_unique_s2 <- geneAssign_df_annot_func %>%
dplyr::filter(grepl("0-.-.-1-.", Strain_interaction)) %>%
dplyr::select(sseqid, ko_name, ko_function_abbrev:ko_level_C) %>%
dplyr::filter(!is.na(ko_name)) %>%
unique() %>%
ggplot(., aes(x = ko_level_C, fill = ko_level_B))+
# geom_point()+
# geom_boxplot()+
geom_bar(stat = "count", col = "black")+
scale_fill_brewer("", palette = "Paired")+
theme_bw()+
theme(axis.text.x=element_text(size=12, angle = 30, hjust = 1),
axis.text.y = element_text(size = 12),
axis.title=element_text(size=16),
title=element_text(size=16), legend.text=element_text(size=16),
strip.text = element_text(size = 16))+
guides(fill = FALSE)
# ylim(0,40)
cowplot::plot_grid(p_unique_s1, p_unique_s2, align = "hv", ncol = 2)
Gene enrichment
# Test for enrichment of functional categories in strains with small inferred genomes
## Test for enrichment of KO level C terms
bg_gsea <- geneAssign_df_annot_func %>%
dplyr::select(ko_level_C, sseqid) %>%
distinct()
bg_gsea[is.na(bg_gsea)] <- "Unknown"
s1g <- geneAssign_df_annot_func %>%
dplyr::filter(`Strain 5` == 1) %>%
dplyr::select(sseqid) %>%
unique()
strain1_gsea <- enricher(gene = s1g$sseqid,
universe = bg_gsea$sseqid,
TERM2GENE = bg_gsea,
pvalueCutoff = 0.05,
qvalueCutoff = 0.2)Strain DOM usage
# Import DOM usage table
DOM_usage_df <- read.csv("./IMG_annotation/DOM_usage.csv", stringsAsFactors = FALSE,
check.names =FALSE)
DOM_usage_df$DOM_type <- gsub("\\n","\n", DOM_usage_df$DOM_type, fixed = TRUE)
# Retain COG_ids that are found in DOM_usage_df list
COG_profiles_sub <- geneAssign_df_annot_func %>%
dplyr::filter(cog_id %in% DOM_usage_df$COG_ID)
COG_profiles_sub <- dplyr::left_join(COG_profiles_sub, DOM_usage_df,
by = c("cog_id" = "COG_ID")) %>%
dplyr::select(sseqid, Strain_interaction, cog_id, DOM_type) %>%
unique()
# Add COG counts
COG_profiles_sub_cnt <- COG_profiles_sub %>%
group_by(Strain_interaction, DOM_type) %>%
mutate(sum_counts = n())
# Order COG ids according to classification
COG_order <- DOM_usage_df %>%
dplyr::filter(COG_ID %in% unique(COG_profiles_sub_cnt$cog_id))
COG_profiles_sub_cnt$cog_id <- factor(COG_profiles_sub_cnt$cog_id,
levels = COG_order$COG_ID)
# # make heatmap
# hm_DOC1 <- complete(COG_profiles_sub_cnt, cog_id,
# fill = list(sum_counts = 0)) %>%
# # dplyr::filter(Genome %in% selected_genomes) %>%
# ggplot(aes(y = Strain_interaction, x= cog_id)) + # x and y axes => Var1 and Var2
# geom_tile(aes(fill = sum_counts), col = "lightgrey") + # background colours are mapped according to the value column
# geom_text(aes(label = round(sum_counts, 0)), size = 3) + # write the values
# # scale_fill_gradientn(colours = terrain.colors(10), trans = "log1p")+
# # scale_fill_gradient(low = "lightblue", high = "darkslategray", na.value="white",
# # trans = "log1p", limits=c(1, 40)) +
# scale_fill_distiller(palette="YlOrRd", na.value="lightgrey", trans = "sqrt",
# direction = 1) +
# theme(panel.grid.major.x=element_blank(), #no gridlines
# panel.grid.minor.x=element_blank(),
# panel.grid.major.y=element_blank(),
# panel.grid.minor.y=element_blank(),
# panel.background=element_rect(fill="white"), # background=white
# axis.text.x = element_text(angle=45, hjust = 1, vjust=1, size = 12,face = "bold"),
# plot.title = element_text(size=20,face="bold"),
# axis.text.y = element_text(size = 12,face = "bold"))+
# theme(legend.title=element_text(face="bold", size=14)) +
# scale_x_discrete(name="") +
# scale_y_discrete(name="") +
# labs(fill="Gene\ncount")
# print(hm_DOC1)
# Barplot
p_DOM_dist <- COG_profiles_sub_cnt %>%
dplyr::filter(Strain_interaction != "0-0-0-0-0") %>%
dplyr::select(Strain_interaction, DOM_type, sum_counts) %>%
unique() %>%
ggplot2::ggplot(aes(x = Strain_interaction, y = sum_counts, fill = DOM_type))+
geom_bar(alpha = 0.8, stat = "identity", color = "black", width = 0.5)+
scale_fill_brewer("", palette = "Paired")+
theme_bw()+
theme(axis.text=element_text(size=14), axis.title=element_text(size=20),
title=element_text(size=20), legend.text=element_text(size=10),
legend.background = element_rect(fill="transparent"),
# axis.text.x = element_text(angle = 65, hjust = 1),
strip.text.x=element_text(size = 18),
legend.position="right",
axis.text.x=element_text(size = 14, angle =45, hjust= 1),
plot.title = element_text(hjust = 0, size=18))+
guides(fill=guide_legend(ncol=1))+
ylab("Number of genes")+
xlab('Intersection')
print(p_DOM_dist)
9. iRep analysis
# Import data
irep_files <- list.files("./iRep/results", pattern = "irep_")
irep_sample_order <- read.table("./iRep/results/sample_order.tsv")
for(file in irep_files){
# print(file)
file_tmp <- read.table(paste("./iRep/results/", file, sep = ""),
stringsAsFactors = FALSE)
bin_name <- gsub("irep_results_|.tsv", "", file)
# print(bin_name)
file_tmp <- data.frame(Genome_bin = bin_name, iRep_value = file_tmp[,2],
iRep_type = rep(c("iRep_filtered","iRep_unfiltered"), 24),
Sample = rep(irep_sample_order$V1, each = 2))
if(file == irep_files[1]) results_irep <- file_tmp else{
results_irep <- rbind(results_irep, file_tmp)
}
}
# Rename sample names
results_irep$Sample <- gsub(".A", "", results_irep$Sample, fixed = TRUE)
results_irep$Sample <- gsub(".C", "", results_irep$Sample, fixed = TRUE)
# Rename genome bins
results_irep <- dplyr::left_join(results_irep, new_bin_names, by = c("Genome_bin" = "bins"))
results_irep$new_bin_name <- as.character(results_irep$new_bin_name)
results_irep$new_bin_name <- factor(results_irep$new_bin_name, levels =
c("MAG1.FA-MLB-DN","MAG2.FA-MLB-SN",
"MAG3.FA-MLB-SN", "MAG4.FA-M110-DN",
"MAG5.SP-M110-DD","MAG6.SP-M15-SD",
"MAG7.SU-MLB-SD","MAG8.SU-M110-DCMD",
"MAG9.SU-M15-SN","MAG10.SU-M15-SN"))
# Merge with metadata
results_irep <- dplyr::left_join(results_irep, meta, by = "Sample")
results_irep$iRep_value <- as.numeric(as.character(results_irep$iRep_value))
# Format metadata
results_irep$Site <- as.character(results_irep$Site)
results_irep$Site <- gsub("110", "Lake Michigan\nsite M110", results_irep$Site)
results_irep$Site <- gsub("15", "Lake Michigan\nsite M15", results_irep$Site)
results_irep$Site <- gsub("Buoy", "Muskegon Lake", results_irep$Site)
results_irep$Site <- factor(results_irep$Site, levels = c("Muskegon Lake",
"Lake Michigan\nsite M15",
"Lake Michigan\nsite M110"))
results_irep$Season <- as.character(results_irep$Season)
results_irep$Season <- factor(results_irep$Season, levels = c("Spring", "Summer","Fall"))Plot results
# Plot unfiltered iRep
irep_p1 <- results_irep %>%
dplyr::filter(iRep_type == "iRep_unfiltered") %>%
ggplot(aes(x = new_bin_name, y = iRep_value, fill = new_bin_name, col = new_bin_name))+
geom_point(color = "black", alpha = 0.7, size = 4, shape = 21)+
# scale_shape_manual(values = c(21,22,23))+
geom_boxplot(alpha = 0.3, col = "black", outlier.shape = NA, size = 0.3)+
scale_fill_manual("", values = fill_palette)+
scale_color_brewer("",palette = "Accent")+
theme_bw()+
# geom_point(size = 4, color = "black", alpha = 0.7)+
# scale_shape_manual(values = c(21,24,23))+
# geom_boxplot(alpha = 0.4)+
theme(axis.text=element_text(size=12), axis.title=element_text(size=16),
title=element_text(size=16), legend.text=element_text(size=10),
legend.background = element_rect(fill="transparent"),
axis.text.x = element_text(size = 12, angle = 45, hjust = 1),
strip.text=element_text(size=12), legend.position = "bottom",
strip.background = element_rect(fill = adjustcolor("gray", 0.15)))+
ylab(paste0("iRep"))+
guides(fill=FALSE)+
facet_grid(Season~Site, scales ="free")+
xlab("")+
scale_y_continuous(labels=scaleFUN, limits = c(1,4.5))+
# geom_smooth(se = FALSE)+
guides(fill=FALSE)
# coord_trans(y = "sqrt")
print(irep_p1)
# Plot filtered irep
irep_p2 <- results_irep %>%
dplyr::filter(iRep_type == "iRep_filtered") %>%
ggplot(aes(x = new_bin_name, y = iRep_value, fill = new_bin_name, col = new_bin_name))+
geom_point(color = "black", alpha = 0.7, size = 4, shape = 21)+
# scale_shape_manual(values = c(21,22,23))+
geom_boxplot(alpha = 0.3, col = "black", outlier.shape = NA, size = 0.3)+
scale_fill_manual("", values = fill_palette)+
scale_color_brewer("",palette = "Accent")+
theme_bw()+
# geom_point(size = 4, color = "black", alpha = 0.7)+
# scale_shape_manual(values = c(21,24,23))+
# geom_boxplot(alpha = 0.4)+
theme(axis.text=element_text(size=12), axis.title=element_text(size=16),
title=element_text(size=16), legend.text=element_text(size=10),
legend.background = element_rect(fill="transparent"),
axis.text.x = element_text(size = 12, angle = 45, hjust = 1),
strip.text=element_text(size=12), legend.position = "bottom",
strip.background = element_rect(fill = adjustcolor("gray", 0.15)))+
ylab(paste0("iRep"))+
guides(fill=FALSE)+
facet_grid(Season~Site, scales ="free")+
xlab("")+
scale_y_continuous(labels=scaleFUN, limits = c(1,4.5))+
# geom_smooth(se = FALSE)+
guides(fill=FALSE)
# coord_trans(y = "sqrt")
print(irep_p2)
# Only focus on MAGs for which we can determine accurate iRep values
irep_p3 <- results_irep %>%
dplyr::filter(iRep_type == "iRep_filtered") %>%
dplyr::filter(!is.na(iRep_value)) %>%
ggplot(aes(x = Season, y = iRep_value, fill = new_bin_name, col = new_bin_name))+
geom_point(color = "black", alpha = 0.7, size = 4, shape = 21)+
# scale_shape_manual(values = c(21,22,23))+
geom_boxplot(alpha = 0.3, col = "black", outlier.shape = NA, size = 0.3)+
scale_fill_manual("", values = fill_palette)+
scale_color_brewer("",palette = "Accent")+
theme_bw()+
# geom_point(size = 4, color = "black", alpha = 0.7)+
# scale_shape_manual(values = c(21,24,23))+
# geom_boxplot(alpha = 0.4)+
theme(axis.text=element_text(size=12), axis.title=element_text(size=16),
title=element_text(size=16), legend.text=element_text(size=10),
legend.background = element_rect(fill="transparent"),
axis.text.x = element_text(size = 12, angle = 45, hjust = 1),
strip.text=element_text(size=12), legend.position = "bottom",
strip.background = element_rect(fill = adjustcolor("gray", 0.15)))+
ylab(paste0("iRep"))+
guides(fill=FALSE)+
facet_grid(new_bin_name~Site, scales ="free")+
xlab("")+
scale_y_continuous(labels=scaleFUN, limits = c(1,4.5))+
# geom_smooth(se = FALSE)+
guides(fill=FALSE)
# coord_trans(y = "sqrt")
print(irep_p3)
# Only focus on MAGs for which we can determine accurate iRep values
irep_p4 <- results_irep %>%
dplyr::filter(iRep_type == "iRep_filtered") %>%
dplyr::filter(!is.na(iRep_value)) %>%
ggplot(aes(x = new_bin_name, y = iRep_value, fill = Site, col = Site))+
# geom_point(color = "black", alpha = 0.7, size = 4, shape = 21)+
# scale_shape_manual(values = c(21,22,23))+
geom_boxplot(alpha = 0.3, col = "black", outlier.shape = NA, size = 0.3)+
scale_fill_brewer("", palette = "Accent")+
scale_color_brewer("",palette = "Accent")+
theme_bw()+
# geom_point(size = 4, color = "black", alpha = 0.7)+
# scale_shape_manual(values = c(21,24,23))+
# geom_boxplot(alpha = 0.4)+
theme(axis.text=element_text(size=12), axis.title=element_text(size=16),
title=element_text(size=16), legend.text=element_text(size=10),
legend.background = element_rect(fill="transparent"),
axis.text.x = element_text(size = 12, angle = 45, hjust = 1),
strip.text=element_text(size=12), legend.position = "bottom",
strip.background = element_rect(fill = adjustcolor("gray", 0.15)))+
ylab(paste0("iRep"))+
# facet_grid(new_bin_name~Site, scales ="free")+
xlab("")+
scale_y_continuous(labels=scaleFUN, limits = c(1,4.5))
# geom_smooth(se = FALSE)+
# guides(fill=FALSE)
# coord_trans(y = "sqrt")
print(irep_p4)
# Focus on all MAGs
irep_p5 <- results_irep %>%
dplyr::filter(iRep_type == "iRep_unfiltered") %>%
dplyr::filter(!is.na(iRep_value)) %>%
ggplot(aes(x = new_bin_name, y = iRep_value, fill = Site, col = Site))+
# geom_point(color = "black", alpha = 0.7, size = 4, shape = 21)+
# scale_shape_manual(values = c(21,22,23))+
geom_boxplot(alpha = 0.3, col = "black", outlier.shape = NA, size = 0.3)+
scale_fill_brewer("", palette = "Accent")+
scale_color_brewer("",palette = "Accent")+
theme_bw()+
# geom_point(size = 4, color = "black", alpha = 0.7)+
# scale_shape_manual(values = c(21,24,23))+
# geom_boxplot(alpha = 0.4)+
theme(axis.text=element_text(size=12), axis.title=element_text(size=16),
title=element_text(size=16), legend.text=element_text(size=10),
legend.background = element_rect(fill="transparent"),
axis.text.x = element_text(size = 12, angle = 45, hjust = 1),
strip.text=element_text(size=12), legend.position = "bottom",
strip.background = element_rect(fill = adjustcolor("gray", 0.15)))+
ylab(paste0("iRep"))+
# facet_grid(new_bin_name~Site, scales ="free")+
xlab("")+
scale_y_continuous(labels=scaleFUN, limits = c(1,4.5))
# geom_smooth(se = FALSE)+
# guides(fill=FALSE)
# coord_trans(y = "sqrt")
print(irep_p5)
Abundances
# Merge irep results with abundance estimates
# Correlate abundances with iRep estimates
# PlotPangenome analysis
# Import and format panG data
# Need to use fread because of the way additional COG categories are interspersed by tabs...
panG <- data.table::fread("./panG/SUMMARY_Final_bins_wMIXED/panG-LIM-ALL_protein_clusters_summary.txt")
# Remove COG annotation
panG <- panG[, -c("V1","V2","V6", "V7")]
# Add column names
colnames(panG) <- c("bin_name", "genome_name", "gene_callers_id",
"COG_FUNCTION_ACC", "COG_FUNCTION", "aa_sequence")
panG_MAG <- panG
panG_MAG$gene_callers_id <- as.character(panG_MAG$gene_callers_id)
# Add unique gene identifier
panG_MAG$unique_gene_callers_id <- interaction(panG_MAG$genome_name,
panG_MAG$gene_callers_id,
sep = ".")
sum(panG_MAG$genome_name=="MAG8_SU_M110_DCMD")## [1] 2294panG_MAG$aa_sequence <- gsub("-", "", panG_MAG$aa_sequence)
# Export all AA sequences for annotation with KAAS or for blast
# for(i in 1:length(unique(panG_MAG$genome_name))){
# tmp <- panG_MAG %>% dplyr::filter(genome_name == unique(panG_MAG$genome_name)[i]) %>%
# droplevels()
# write.table(file = paste0("./panG/", unique(panG_MAG$genome_name)[i],
# "_aa_export.fa"),
# paste0(">", tmp$unique_gene_callers_id,
# "\n", tmp$aa_sequence, sep = ""),
# quote = FALSE, row.names = FALSE, col.names = FALSE
# )
# }
# Import KEGG annotation through KAAS (http://www.genome.jp/tools/kaas/) of amino
# acid sequences
ko_files <- list.files(".", pattern = "KO_annot",
recursive = TRUE)
panG_ko <- data.frame()
for(ko_file in ko_files){
tmp <- data.table::fread(ko_file, fill = TRUE, sep = "\t", header = FALSE)
tmp <- data.frame(tmp, Genome = do.call(rbind, strsplit(ko_file, "-"))[,2])
if(ko_file == ko_files[1]) panG_ko <- tmp else{
panG_ko <- rbind(panG_ko, tmp)
}
}
panG_ko$Genome <- gsub(".txt","",panG_ko$Genome)
colnames(panG_ko)[1:2] <- c("unique_gene_callers_id", "ko_id")
panG_ko <- panG_ko[panG_ko$ko_id != "",]
panG_ko$ko_id <- gsub(" ","", panG_ko$ko_id)
# Annotate KO_IDs with KEGG hierarchy
panG_ko <- dplyr::left_join(panG_ko, ko_path_df, by = "ko_id")
# join with corresponding COG ids
panG_ko_cog <- dplyr::left_join(panG_MAG, panG_ko,
by = "unique_gene_callers_id")
# Shorten/change ko_level_B annotation a bit
panG_ko_cog$ko_level_B[panG_ko_cog$ko_level_B == "Cellular community - prokaryotes"] <- "Biofilm formation & quorum sensing"
panG_ko_cog$ko_level_B[panG_ko_cog$ko_level_B == "Xenobiotics biodegradation and metabolism"] <- "Xenobiotics degradation"
panG_ko_cog$ko_level_C[panG_ko_cog$ko_level_C == "Biofilm formation - Escherichia coli "] <- "Biofilm formation"
panG_ko_cog$ko_level_C[panG_ko_cog$ko_level_C == "Biofilm formation - Pseudomonas aeruginosa "] <- "Biofilm formation"
panG_ko$ko_level_C[panG_ko$ko_level_C == "Biofilm formation - Escherichia coli "] <- "Biofilm formation"
panG_ko$ko_level_C[panG_ko$ko_level_C == "Biofilm formation - Pseudomonas aeruginosa "] <- "Biofilm formation"
# Add column denoting whether it is core/mixed or accessory
panG_ko_cog$bin_core <- factor(panG_ko_cog$bin_name == "CORE_PC" | panG_ko_cog$bin_name == "MIXED_PCs")
panG_ko_cog$bin_core <- plyr::revalue(panG_ko_cog$bin_core, replace = c("TRUE" = "CORE/Mixed", "FALSE" = "Accessory"))# Perform hypergeometric test to see if there is functional enrichment in the pangenome
for(i_genome in 1:length(unique(panG_ko$Genome))){
bg_gsea <- panG_ko %>% dplyr::filter(Genome == unique(panG_ko$Genome)[i_genome]) %>%
dplyr::select(ko_level_C,unique_gene_callers_id) %>%
distinct()
panG_gseq <- panG_ko_cog %>%
dplyr::filter(Genome == unique(panG_ko$Genome)[i_genome] &
bin_core == "Accessory") %>%
dplyr::select(ko_level_C,unique_gene_callers_id) %>%
distinct()
panG_gsea <- enricher(gene = panG_gseq$unique_gene_callers_id,
universe = bg_gsea$unique_gene_callers_id,
TERM2GENE = bg_gsea,
pvalueCutoff = 0.05,
qvalueCutoff = 0.2)
tmp_gseq <- data.frame(panG_gsea@result, Genome = unique(panG_ko$Genome)[i_genome])
if(i_genome == 1) results_gseq <- tmp_gseq else {
results_gseq <- rbind(results_gseq, tmp_gseq)
}
genes_per_gseq <- panG_gseq %>%
dplyr::filter(ko_level_C %in% panG_gsea@result$Description) %>%
group_by(ko_level_C) %>%
summarize(Counts = n())
genes_per_gseq <- data.frame(genes_per_gseq,
Genome = unique(panG_ko$Genome)[i_genome])
if(i_genome == 1) results_pergene <- genes_per_gseq else {
results_pergene <- rbind(results_pergene, genes_per_gseq)
}
}
# Calculate fraction of enriched genes vs. annotated panG genes
results_pergene_sum <- results_pergene %>%
group_by(Genome) %>%
summarize(enrich_counts = sum(Counts))
annotated_fraction_panG <- panG_ko_cog %>%
dplyr::filter(bin_core == "Accessory") %>%
group_by(genome_name) %>%
dplyr::select(unique_gene_callers_id) %>%
distinct() %>%
summarize(panG_counts = n())## Adding missing grouping variables: `genome_name`print(annotated_fraction_panG)## # A tibble: 19 x 2
## genome_name panG_counts
## <chr> <int>
## 1 MAG1_FA_MLB_DN 551
## 2 MAG10_SU_M15_SN 658
## 3 MAG2_FA_MLB_SN 1020
## 4 MAG3_FA_MLB_SN 537
## 5 MAG4_FA_M110_DN 647
## 6 MAG5_SP_M110_DD 277
## 7 MAG6_SP_M15_SD 554
## 8 MAG7_SU_MLB_SD 120
## 9 MAG8_SU_M110_DCMD 288
## 10 MAG9_SU_M15_SN 234
## 11 REF_Lim_103DPR2 333
## 12 REF_Lim_2KL27 192
## 13 REF_Lim_2KL3 674
## 14 REF_Lim_63ED37_2 200
## 15 REF_Lim_DM1 709
## 16 REF_Lim_IID5 560
## 17 REF_Lim_Rim11 287
## 18 REF_Lim_Rim28 447
## 19 REF_Lim_Rim47 187results_pergene_sum <- left_join(results_pergene_sum, annotated_fraction_panG,
by = c("Genome" = "genome_name"))
# Merge actual unique gene counts per gsea category with the accessory genome sizes
results_pergene_sum <- results_pergene_sum %>%
mutate(frac_enrich = enrich_counts/panG_counts)
print(results_pergene_sum)## # A tibble: 10 x 4
## Genome enrich_counts panG_counts frac_enrich
## <chr> <int> <int> <dbl>
## 1 MAG1_FA_MLB_DN 87 551 0.158
## 2 MAG10_SU_M15_SN 91 658 0.138
## 3 MAG2_FA_MLB_SN 180 1020 0.176
## 4 MAG3_FA_MLB_SN 53 537 0.0987
## 5 MAG4_FA_M110_DN 144 647 0.223
## 6 MAG5_SP_M110_DD 14 277 0.0505
## 7 MAG6_SP_M15_SD 32 554 0.0578
## 8 MAG7_SU_MLB_SD 15 120 0.125
## 9 MAG8_SU_M110_DCMD 66 288 0.229
## 10 MAG9_SU_M15_SN 5 234 0.0214# Make plot
results_gseq$Genome <- gsub("_","-",results_gseq$Genome)
results_gseq$Genome <- gsub("-FA",".FA",results_gseq$Genome)
results_gseq$Genome <- gsub("-SU",".SU",results_gseq$Genome)
results_gseq$Genome <- gsub("-SP",".SP",results_gseq$Genome)
results_gseq$Genome <- factor(results_gseq$Genome, levels =
c("MAG5.SP-M110-DD","MAG2.FA-MLB-SN",
"MAG3.FA-MLB-SN", "MAG4.FA-M110-DN",
"MAG1.FA-MLB-DN", "MAG10.SU-M15-SN",
"MAG6.SP-M15-SD", "MAG8.SU-M110-DCMD",
"MAG7.SU-MLB-SD","MAG9.SU-M15-SN"))
p_panG5 <- results_gseq %>%
ggplot(aes(x = Description, fill = Genome, y = Count))+
geom_bar(color = "black", stat = "identity")+
theme_bw()+
scale_fill_brewer(palette="Paired")+
ylab("Number of genes") + xlab("")+
facet_grid(.~Genome, scales = "free")+
theme(axis.text.y=element_text(size=12.5), axis.title.y =element_text(size = 18),
legend.background = element_rect(fill="transparent"),
axis.text.x = element_text(size=12.5, angle = 45, hjust = 1),
plot.margin = unit(c(1.1,1.1,1.1,1.1), "cm"),
strip.text = element_text(size = 12))+
guides(fill = FALSE)
# print(p_panG5)
# Make heatmap instread
results_gseq_wide <- complete(results_gseq, Genome, Description,
fill = list(Count = 0)) %>%
dplyr::select(Description, Genome, Count) %>%
spread(., Description, Count)
dat <- results_gseq_wide[,2:ncol(results_gseq_wide)] # numerical columns
rownames(dat) <- results_gseq_wide$Genome
# row.order <- hclust(dist(dat))$order # clustering
col.order <- hclust(dist(t(dat)))$order
dat_new <- dat[, col.order] # re-order matrix accoring to clustering
rownames(dat_new) <- results_gseq_wide$Genome
df_molten_dat <- melt(as.matrix(dat_new)) # reshape into dataframe
names(df_molten_dat)[c(1:2)] <- c("Genome", "Description")
p_panG6 <- complete(df_molten_dat, Genome, Description,
fill = list(value = 0)) %>%
# dplyr::filter(Genome %in% selected_genomes) %>%
ggplot(aes(x = Genome, y = Description)) + # x and y axes => Var1 and Var2
geom_tile(aes(fill = value), col = "lightgrey") + # background colours are mapped according to the value column
geom_text(aes(label = round(value, 0)), size = 4) + # write the values
# scale_fill_gradientn(colours = terrain.colors(10), trans = "log1p")+
# scale_fill_gradient(low = "lightblue", high = "darkslategray", na.value="white",
# trans = "log1p", limits=c(1, 40)) +
scale_fill_distiller(palette="YlOrRd", na.value="lightgrey", trans = "sqrt",
direction = 1) +
theme(panel.grid.major.x=element_blank(), #no gridlines
panel.grid.minor.x=element_blank(),
panel.grid.major.y=element_blank(),
panel.grid.minor.y=element_blank(),
panel.background=element_rect(fill="white"), # background=white
axis.text.x = element_text(angle=45,
hjust = 0, vjust=0, size = 12,face = "bold"),
plot.title = element_text(size=20,face="bold"),
axis.text.y = element_text(size = 12,face = "bold"))+
theme(legend.title=element_text(face="bold", size=14)) +
scale_x_discrete(name="", position = "top") +
scale_y_discrete(name="") +
labs(fill="Gene\ncount")
print(p_panG6)
Phenotypic diversity
Plot phenoD
# Plot results
p_MAG_Pdiv2 <- ggplot(results_pd, aes(x = Season, y = D2))+
theme_bw()+
scale_fill_brewer(palette = "Accent")+
geom_point(size = 4, color = "black",
alpha = 0.7, aes(fill = Season, shape = Depth))+
scale_shape_manual(values = c(21,24,23))+
geom_boxplot(alpha = 0.4, width = 0.2, outlier.shape = NA,
aes(fill = Season))+
theme(axis.text=element_text(size=14), axis.title=element_text(size=20),
title=element_text(size=20), legend.text=element_text(size=12),
legend.background = element_rect(fill="transparent"),
axis.text.x = element_text(size = 14, angle = 45, hjust = 1),
strip.text=element_text(size=14), legend.position = "bottom",
strip.background = element_rect(fill = adjustcolor("gray", 0.15)))+
ylab(expression("Transcriptional diversity (D"[2]*")"))+
guides(
fill = guide_legend(override.aes=list(colour=brewer.pal(3,"Accent"))))+
facet_grid(Site~new_bin_name, scales ="free")+
xlab("")+
# geom_errorbar(aes(ymin = D2 - sd.D2, ymax = D2 + sd.D2), width = 0.05)+
scale_y_continuous(labels=scaleFUN)
# coord_trans(y = "sqrt")
# print(p_MAG_Pdiv2)
# Plot results
dodge <- position_dodge(width=0.6)
p_MAG_Pdiv3 <- ggplot(results_pd, aes(x = new_bin_name, y = D2))+
theme_bw()+
scale_fill_brewer("",palette = "Accent")+
# geom_point(size = 4, color = "black",
# alpha = 0.5, aes(fill = Season, shape = Depth),
# position=dodge)+
geom_boxplot(alpha = 0.4, width = 0.4, outlier.shape = NA,
aes(fill = Season),
position=dodge)+
scale_shape_manual("",values = c(21,24,23))+
# geom_boxplot(alpha = 0.4, width = 0.2, outlier.shape = NA)+
theme(axis.text=element_text(size=14), axis.title=element_text(size=20),
title=element_text(size=20), legend.text=element_text(size=12),
legend.background = element_rect(fill="transparent"),
axis.text.x = element_text(size = 14, angle = 45, hjust = 1),
strip.text=element_text(size=14), legend.position = "bottom",
strip.background = element_rect(fill = adjustcolor("gray", 0.15)))+
ylab(expression("Transcriptional diversity (D"[2]*")"))+
guides(
fill = guide_legend(override.aes=list(colour=brewer.pal(3,"Accent"))))+
# facet_grid(Season~., scales ="free")+
xlab("")+
# geom_errorbar(aes(ymin = D2 - sd.D2, ymax = D2 + sd.D2), width = 0.05)+
scale_y_continuous(labels=scaleFUN)
# coord_trans(y = "sqrt")
print(p_MAG_Pdiv3)
Match with abundance
# Merge with abundance data
sel_abund <- c("new_bin_name","Sample",
"norm_lower_rel_abundance",
"norm_mean_rel_abundance",
"norm_upper_rel_abundance",
"MGT","Topt","Lineage")
results_pd_abund <- dplyr::left_join(results_pd, data_total_abund[, sel_abund],
by = c("new_bin_name", "Sample"))
p_MAG_Pdiv4 <- ggplot(results_pd_abund, aes(y = D2, x = norm_mean_rel_abundance))+
theme_bw()+
geom_point(size = 4, color = "black", shape = 21,
alpha = 0.5, aes(fill = new_bin_name))+
scale_fill_brewer("",palette = "Paired")+
# geom_boxplot(alpha = 0.4, width = 0.2, outlier.shape = NA)+
theme(axis.text=element_text(size=14), axis.title=element_text(size=20),
title=element_text(size=20), legend.text=element_text(size=12),
legend.background = element_rect(fill="transparent"),
axis.text.x = element_text(size = 14),
strip.text=element_text(size=14), legend.position = "bottom",
strip.background = element_rect(fill = adjustcolor("gray", 0.15)))+
ylab(expression("Transcriptional diversity (D"[2]*")"))+
xlab("Mean normalized abundance")+
# guides(
# fill = guide_legend(override.aes=list(colour=brewer.pal(10,"Paired"))))
# facet_grid(Season~Site, scales ="free")
# geom_errorbar(aes(ymin = D2 - sd.D2, ymax = D2 + sd.D2), width = 0.05)+
coord_trans(y = "log10", x = "log10")
print(p_MAG_Pdiv4)
p_MAG_Pdiv5 <- ggplot(results_pd_abund, aes(y = D2, x = log(2)/MGT))+
theme_bw()+
geom_point(size = 4, color = "black", shape = 21,
alpha = 0.5, aes(fill = new_bin_name))+
scale_fill_brewer("",palette = "Paired")+
# geom_boxplot(alpha = 0.4, outlier.shape = NA, aes(fill = new_bin_name),
# width = 0.005)+
theme(axis.text=element_text(size=14), axis.title=element_text(size=20),
title=element_text(size=20), legend.text=element_text(size=12),
legend.background = element_rect(fill="transparent"),
axis.text.x = element_text(size = 14),
strip.text=element_text(size=14), legend.position = "bottom",
strip.background = element_rect(fill = adjustcolor("gray", 0.15)))+
ylab(expression("Transcriptional diversity (D"[2]*")"))+
xlab(expression("Maximum growth rate (µ"^-1*")"))+
stat_summary(fun.data=mean_sdl, fun.args = list(mult = 1),
geom="pointrange", color="#333333", size = 1)+
guides(
fill = guide_legend(override.aes=list(colour=brewer.pal(10,"Paired"))))
# facet_grid(.~Lineage, scales ="free")
# geom_errorbar(aes(ymin = D2 - sd.D2, ymax = D2 + sd.D2), width = 0.05)+
# coord_trans(x = "log10")
print(p_MAG_Pdiv5)